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ABSaUCT 

Tht  problM  of  tht  nuaorloal  solution  by  the  method  of 
eharftcterlsties  of  the  two-  and  three-dimensional  flow  of  an 
Inwlsold  nonequlllbrlua  gas  Is  formulated.  For  unsteady  flows 
the  method  Is  limited^  for  practical  reasons,  to  the  solution 
of  hypersonic  or  other  flows  which  have  short  transient  tines. 
The  specific  details  of  the  finite  difference  network  to  be 
used  and  the  associated  finite  difference  equations  are 
presented.  It  was  found  that  numerical  Instability  was 
possible  with  certain  previously  proposed  networks.  An 
existing  stability  criterion  was  found  to  apply  to  the  multi¬ 
dimensional  method  of  characteristics  and  this  criterion  wae 
used  to  synthesise  a  stable  finite'  difference  network. 

Ugh  speed  digital  oosputer  programs  are  presented  which 
perform  the  numerical  calculation.  Because  the  operating 
speed  and  storage  capacity  of  the  larger  and  faster  computers 
available  today  does  not  allow  the  solution  of  the  most 
general  problem  In  a  reasonable  length  of  time.  It  was 
necessary  to  limit  the  programmed  calculation  to  the  two- 
dimensional  unsteady  flow  of  a  perfect  gas.  Rather  long 
programs  which  require  large  amounts  of  Input  data  and  produce 
even  larger  amounts  of  output  data  are  the  results  of  the 
programming  effort. 

Four  exaiqple  oases  of  the  flow  about  a  circular  cylinder 
with  Its  axis  perpendicular  to  a  Mach  five  free  stream  are 
presented.  The  cylinder  is  held  steady,  symmetrically  and 
asypmetrlcally  warped  to  an  elliptic  cross-section  and 
oscillated  in  a  direction  normal  to  the  free  stream.  These 
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oattf  Indicate  that  the  aolutlon  of  ■ultl-dlaenalonal  flowa 
by  the  aethod  of  oharaeterlatloa  la  both  practical  and 
feaalble.  With  the  aubatltutlon  of  higher  order  Intexpolatlon 
for  the  linear  Intexpolatlon  used  In  thla  Initial  atndy, 
accurate  and  uaeful  reaulta  for  engineering  probleaa  can 
be  obtained. 

The  extension  of  the  progranned  procedures  to  the  most 
general  problem  which  can  be  solved  with  the  method  of 
characteristics  Is  discussed  and  the  steps  required  to 
extend  the  procedxire  are  Indicated  • 
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CHAPTER  1 


INTRODUCTION 


This  study  was  Inltlatad  with  tht  purpose  of  analysing 
inviscid  hypersonic  flow  about  three-diaensional  bodies  in 
unsteady  motion*  In  the  course  of  this  investigation,  a  more 
general  method  has  been  considered  which  applies  to  all  inviscid 
unsteady  flows,  not  Just  hypersonic  flows*  There  are  practical 
limitations  on  the  application  of  the  method,  however* 

The  study  of  hypersonic  flow  is  mathematically  nonlinear 
in  most  cases,  perhaps  by  definition*  A  review  of  the  litera¬ 
ture  reveals  that  because  of  this  nonlinearity,  theoretical  study 
is  very  difficult*  Many  approximations  are  usually  made  so  that 
the  mathematical  problem  becomes  amenable  to  analytical  methods* 
See,  for  example,  the  book  by  Hayes  and  Probstein^*  Many  times 
it  is  very  difficult  to  evaluate  the  effect  of  the  approxima¬ 
tions  on  the  results  of  the  analysis  after  the  approximate 
solution  has  been  obtained.  When  unsteady  hypersonic  flow  is 
considered,  the  problems  are  even  more  challenging.  Here,  by 
unsteady  hypersonic  flow,  we  mean  truly  unsteady  flow  where  the 
"quasi- steady"  approximation  which  simplifies  certain  unsteady 
hjrpersonic  analysis  is  not  valid.  Also,  without  exact  solutions 
against  which  these  approximate  theories  can  be  checked,  the 


^  Superscripts  refer  to  the  reference  numbers  listed  in  the 
section  entitled  References. 
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•rrori  introduced  by  thea  ere  very  difficult,  if  not  impossible 
to  determine*  Hence,  e  meens  of  exactly  calculating  unsteady 
hypersonic  flow  is  very  desirable* 

The  question  which  faces  us  then  ii,  what  is  the  best  method 
for  determining  exact  solutions  to  unsteady  hypersonic  flows?  A 
method  is  sought  to  solve  the  general  nonlinear  equations  which 
describe  the  flow  of  an  inviscid  nonf^quilibrium  gas*  It  is 
probably  safe  to  assume  that  a  general  analytical  solution  can 
not  be  found  because  of  the  complicated  nature  of  the  equations* 
It  must  be  concluded  then,  that  the  only  hope  lies  in  the  direc¬ 
tion  of  a  numerical  solution* 

# 

The  general,  unsteady  flow  of  a  gas  is  governed  by  a  set  of 
hyperbolic,  quasi- linear  partial  differential  equations*  In 
general,  there  are  two  methods  for  numerically  solving  such 
■ystems  of  equations*  One  method  replaces  thr  derivatives  in 
the  original  set  of  equations  with  finite  differences,  and 
seeks  to  solve  numerically,  this  set  of  difference  equations* 

This  approach,  we  might  call  the  finite  difference  method*  A 
second  procedure  first  writes  the  equations  in  characteristic 
form  to  somewhat  simplify  them  and  then  utilizes  finite  dif¬ 
ferences  to  approximate  the  derivatives*  The  numerical  solu¬ 
tion  is  obtained  by  integration  in  the  characteristic  surfaces, 
hence,  this  approach  is  referred  to  as  the  method  of  character¬ 
istics*  A  discussion  and  examples  of  these  two  basic  approaches 
are  given  in  Chapters  15  and  16  of  the  book  edited  by  Ralston 
•ad 
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The  finite  difference  aethods  night  be  further  subdivided 
Into  stenderd  nethods  end  methods  utilising  ertlflclel  vis¬ 
cosity.  In  the  letter  epproech,  en  extre  ertlflclel  viscosity 
tern  (or  In  sone  ceses,  the  reel  viscosity  ten)  Is  edded  to  the 
equetlons  such  thet  free  bounderles,  e.g.  shock  veves,  ere  euto- 

aetlcelly  detenlned  In  the  solution.  This  epproech  hes  been 

3  4 

studied  by  von  Neumenn  end  Rlchta#yer  »  Lex  ,  end  Lex  end 
Wendroff^  eaong  others.  Thus,  In  generel  there  ere  three  beslc 
numerlcel  epproeches  which  night  be  utilised  In  ettecklng  the 
problem.  Hertree^  hes  proposed  e  method  which  combines  the 
stenderd  finite  difference  epproech  with  the  method  of  cherec- 
terlstlcs  melnly  to  obteln  the  solution  on  en  evenly  speced 
grid  while  still  teklng  edventege  of  the  cherecterlstlcs.  Such 
hybrid  epproeches  usuelly  heve  some  of  the  dlsedventeges  of  the 
orlglnel  methods  end  quite  often  require  more  computetlons  then 
the  orlglnel  methods,  so  thet  they  ere  not  considered  here. 

The  edventeges  end  dlsedventeges  of  eech  of  the  three  beslc 
epproeches  ere  coopered  In  the  following  teble.  Sooe  of  the 
properties  given  In  the  teble  were  obtelned  from  the  books  by 
Fox^  (Chepters  17,  18,  26,  27  end  28),  Forsyth  end  Wesow^  (perts 
1  end  4),  Rlchtmyer^  end  Collets^^  (Chepter  IV,  Pert  5)  while 
some  of  the  estlmetes  ere  the  euthor's. 


•rlton  of  Pottiblt  Mtthodi  of  Solution 


Flnlt*  Dlffcrtnct 

Standard  Plnita  Approach  vlth  Hathod  of 

Dlffaranca  Approach  Artificial  Vtacoaltv  Charactariatica 

I*  Hoar  ara  fraa  boundariaa  handlad? 

apacially  traatad  autOMtically  handlad  apacially  treated 

aa  diacontiauitiaa  vith  artificial  via-  aa  diacontinuitiaa 

coaity  and  by  retaining 
conaa^ativa  propartiaa 
in  the  difference 
aquationa 


2*'  Hov  ara  fixed  boundariaa  handled? 

apacially  treated  apacially  treated »  but  apacially  traatad 

aona  difficulty  can 
ariaa  bacauaa  the  arti¬ 
ficial  viacoaity  in- 
craaaaa  the  order  of 
the  aquationa  and  extra 
boundary  conditiona  ara 
available  which  ara  not 
aaay  to  apacify  in  apse 
caaaa 


3.  Vhat  fornulation  of  the  aquationa  ia  uaually  preferred? 

lularian  or  Lagrangian  in  aoaa  Eularian,  uaually 

Lagrangian  probleaa  utiliaad 


4«  la  a  atability  criterion  required? 

yea  yea  yea,  but  it  ia 

uaually  alaplifiad 
by  knovladae  of  the 
location  or  the 
charactariatic 
aurfacaa 


5*  How  aany  independent  variablaa  can  be  traatad? 

three  apace  diaan-  aora  auitad  to  two  three  apace  diaan- 

aiona  aM  tine  independint  variablaa,  aiona  and  tlaa 

but  can  be  ganeraliaed 


k 


Standard  Flnlta 
Dlffaranca  Approach 


Flnlta  Dlffaranca 
Approach  with 
Artificial  VlacoaltY 


Hathod  of 
Charactarlatlca 


6.  Hov  difficult  la  tha  aachlna  programming  task? 


aoaMifhat  coapllcatad 
baeauaa  boundarlaa 
■ust  ba  traatad 
■pac tally 


almplaat 


■oat  difficult 


7.  Hov  Buch  coaputlng  tlsa  la  raqulrad? 


about  tha  aaaa  aa 

flnlta  dlffaranca 
vlth  artificial 
▼lacoalty 


about  tha  aama  aa 
atandard  flnlta 
dlffaranca  approach 


potentially  tha 
qulckaat  sathod 


8.  How  such  BUichlna  atoraga  la  raqulrad? 


about  tha  aaaa  aa 

flnlta  dlffaranca 
with  artificial 
wlacoalty 


about  tha  aaaa  aa 
atandard  flnlta 
dlffaranca  approach 


allghtly  aora  than 
tha  othar  two 
■athoda 


9*  On  what  aort  of  grid  la  tha  aolutlon  obtalnad? 


avanly  apacad 


avenly  apacad 


character la tic 
grid,  probably 
will  raqulra  inter¬ 
polation  of  tha 
final  raault 


10.  Vhat  la  the  accuracy  of  tha  raaulta? 


Intanadlata 

accuracy 


laaat  accurate  In  that 
tha  free  boundarlaa  are 
not  precisely  located 
and  errors  can  ba  intro¬ 
duced  by  not  properly 
specifying  tha  extra 
boundary,  conditions 
Introduced  by  tha 
artificial  viscosity 


potentially  the 
■ost  accurate, 
also  this  method 
of  solution  most 
closely  follows 
tha  physical  nodal 
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FroB  the  cftlaatcs  and  properties  given  in  the  above  table, 
it  is  possible  to  draw  aose  conclusions  about  the  utility  of 
applying  the  wethods  to  various  types  of  problems.  The  first 
property  listed  in  the  table  shows  a  basic  difference  bet«reen 
the  approach  utilizing  artificial  viscosity  and  the  other  two 
approaches.  This  is  the  major  advantage  of  the  artificial  vis¬ 
cosity  method,  but  it  is  at  the  expense  of  points  2  and  10. 

Hence,  if  the  problem  to  be  solved  involved  a  complicated  inter¬ 
action  of  shock  waves,  for  example,  this  might  be  the  best  method 
to  utilize.  However,  if  only  simple  shock  wave  configurations 
were  to  be  considered,  perhaps  one  of  the  other  two  methods 
might  best  be  utilized.  On  points  1,  2,  3,  and  3  the  standard 
finite  difference  approach  (SFDA)  is  similar  to  the  method  of 
characteristics  (HOC),  %fhile  on  points  6,  8,  and  9  the  SFDA 
seess  to  have  the  advantage,  and  on  points  4,  7,  and  10  the  HOC 
seems  to  be  p'^ef erred. 

The  approach  to  be  adopted  here  must  take  into  account  the 
type  of  problems  that  are  to  be  solved.  In  this  study,  we  are 
interested  generally  in  determining  the  external  hypersonic  flow 
about  bodies  performing  unsteady  motions.  At  least  initially, 
and  in  a  majority  of  the  problems  to  be  considered,  the  bodies 
arc  simple  enough  that  rather  simple  shock  configurations  are 
encountered.  Hence,  the  approach  utilizing  artificial  viscosity 
is  not  required.  The  choice  between  the  SFDA  and  the  HOC  is  not 
as  simple.  As  has  been  indicated,  the  SFDA  has  advantages  on 
three  points  and  the  HOC  has  three,  while  four  can  be  considered 
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equal.  Thus,  the  individual  points  must  be  considered  to  choose 
the  better  method  of  approach.  The  SFDA  is  favored  by  points  6, 

8,  and  9,  trhich  mainly  apply  to  the  ease  of  «rriting  the  program 
and  applying  the  method  to  a  computing  machine.  The  MOC,  on  the 
other  hand,  is  favored  by  points  4,  7,  and  10  which  are  advan¬ 
tages  of  the  finished  working  method.  So,  if  ease  of  %rriting 
the  program  and  initially  applying  the  method  is  most  isq>ortant, 
it  would  appear  that  the  SFDA  is  better,  but  if  extra  initial 
effort  can  be  devoted  to  programning^  a  better  irorking  method 
would  seen  to  result  from  the  HOC.  It  is  felt  that,  because 
once  a  suichine  program  is  mritten  it  can  be  used  over  and  over 
to  solve  many  problems,  it  is  worth  an  extra  initial  effort  to 
%n:ite  a  more  accurate  and  quicker  program.  Hence,  the  MOC  is 
adapted  in  this  study.  Note  that  this  choice  is  based  on  some 
rather  crude  estimates  and  the  validity  of  these  conclusions 
ultimately,  can  be  judged  only  after  'working  programs  have  been 
written  and  are  compared. 

The  MOC  was  first  proposed  by  Massau^^  in  1899  and  is  some¬ 
times  referred  to  as  the  Method  of  Nassau.  The  method  was  adopted 

to  solve  hyperbolic  problems  in  fluid  mechanics,  perhaps  first 

12 

in  1929  by  Prandtl  and  Busemann  .  The  MOC  was  first  used  to 
solve  two-dimensional  irrotational  supersonic  steady  flows  and 
one -dimensional  unsteady  flows  of  a  perfect  gas.  These  calcula¬ 
tions  could  be  made  by  hand  and  in  part  graphically.  Necessarily, 
the  first  applications  of  the  method  had  to  be  limited  to  sim¬ 
plified  problems  in  two  independent  variables,  because  of  the 


7 


iaposslbllity  of  doing  norc  aibltldui  problems  with  hand  cal¬ 
culations.  In  the  1940 *s  and  early  1950* s,  when  electronic 
digital  computers  first  became  available,  the  HOC  was  developed 
extensively  for  more  general  problems  involving  two  independent 
variables.  Ferri^^  and  Meyer^^*^  discuss  this  work  and  give 
an  extensive  list  of  references.  The  rotational  flows  of  real 
gases  became  amenable  to  solution  by  the  NOG  when  computers 
were  used  and  in  some  cases,  were  merely  a  matter  of  the  routine 
application  of  general  computer  programs.  In  the  late  1950 's 
and  early  1960's,  computer  development  had  increased  the  storage 
capacity  and  speed  of  operation  of  the  machines  to  the  point 
that  more  ambitious  problems  could  be  attempted  with  the  MOC. 
Froblems  involving  three  independent  variables  and  nonequilib¬ 
rium  thcrm0d3rnasd.es  are  now. being  attempted  for  the  first  time. 
These  new  advances  in  computer  technology  have  made  it  possible 
to  at  least  consider  the  most  general  t3rpes  of  problems  which 
can  be  solved  with  the  HOC. 

Briefly,  the  MOC  is  used  to  solve  Cauchy  problems  for  a  set 
of  quasi-linear,  h3rperbolic  partial  differential  equations.  Data 
are  given  on  an  initial  surface*  together  with  boundary  con¬ 
ditions  on  the  flow  field,  such  as  solid  body  surfaces  and  shock 
waves.  The  solution  is  then  obtained  on  an  adjacent  surface  by 


*  The  terminology  of  a  problem  involving  three  independent 
variables  is  used  hers  to  emphasise  that  problems  of  more 
than  two  independent  variables  can  be  handled  while  still 
avoidinx  the  mere  general,  but  serfaaps  confusing  termin¬ 
ology  of  four  independent  variables,  e.g.  h3rpersurfaces. 


8 


nuMrlcAl  solution  of  tho  squotlons  written  In  charnct  oris  tie 
form  with  tho  partial  dorlvatlvos  raplacod  by  finite  dlf- 
foroscas.  Tho  solution  on  this  now  surface  Is  then  taken  as 
the  Initial  data  and  the  process  repeated  to  obtain  the  solution 
on  the  next  adjacent  surface.  In  this  way,  the  solution  Is 
obtained  by  marching  In  space  or  time  until  the  desired  solu¬ 
tion  Is  obtained  or  until  the  complete  region  specified  by 
the  Initial  data  has  been  determined. 

A  review  of  the  literature  of  the  HOC  for  two  Independent 
variables  cannot  be  given  here, as  the  volume  of  such  litera¬ 
ture  Is  enormous.  An  attempt  will  be  made  here  to  list  only 
the  most  recent  work  on  the  MCX!  which  involves  problems  which 
Include  more  than  two  Independent  variables  and  those  that  con¬ 
sider  nonequilibrium  thermod3mamlcs. 

The  work  on  problems  with  more  than  two  independent  varia¬ 
bles  can  be  divided  into  two  groups.  One  group,  theoretically 
dtudles  the  problems  and  formulates  then,  and  in  some  cases  even 
goes  so  far  as  to  write  the  finite  difference  equations,  while 
the  second  group  Includes  work  which  has  gone  on  to  actually 
do  some  nusMrlcal  calculations  on  digital  computers.  The  for¬ 
mer  group  Is  by  far  the  larger  of  the  two  groups,  while  the 
latter  has  less  than  ten  members  to  be  sure.  The  reason  for 
this,  is  that  much  time  and  effort  are  required  to  carry  a 
formulation  through  to  a  working  computer  program  or  to  ac¬ 
complish  the  calculation  by  hand  In  simple  cases.  Thus,  many 
studies  of  the  MCX!  have  been  terminated  after  formulating  the 
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problaa.  Included  in  Che  first  group  ere  Che  besic  works  of 

Thomhill^^  end  Cobum  end  Dolph^^.  Seuer^^,  Cllpplnger  end 
19  20 

Giese  end  Holt  heve  elso  suide  conCribuClons  Co  Che  theory. 

21  22 

Cobum  end  Seinl  heve  used  Censor  enelysls  to  fonulece 

e  such  aore  generel  problen  which  hes  Che  perCicuIer  problem 

23 

considered  here  es  e  speciel  cese.  Most  recently,  Fowell 
24 

end  Seuer  heve  published  work  in  this  eree.  In  Che  second 

group  of  pepers,  which  did  sosie  ecCuel  numericel  celculeCions, 

Che  first  few  studies  were  done  when  conputers  were  not  reedily 

eveileble  so  Chet  they  involved  sinplified  celculeCions  which 

could  be  eccoaplished  by  hend  (by  hend  celculeCions,  we  elso 

■een  with  Che  eid  of  desk  celculeCors)  •  These  iniciel  celcule- 

25  26 

Cions  include  Che  work  of  Moeckel  ,  Ferreri  ,  end  Bruhn  end 

27  28  29 

Heeck  •  LeCer  work  by  Butler  end  Tsung  wes  ecconplished 

using  digiCel  cosqimCers.  Results  obCeined  with  the  lecesc 

generecion  of  caBq;>uCers  heve  just  begun  Co  eppeer.  MorreCCi 
30 

et.el.  heve  published  sosm  recent  results.  More  end  isore 
ectuel  celculeCions  ere  being  cerried  out  ec  the  present  else 
on  sore  end  sore  eid>icious  problems  beceuse  of  the  increesed 
verseCiliCy  of  digiCel  cosq>uCers. 

The  study  of  hypersonic  flows  must  ellow  for  Che  pos¬ 
sibility  of  high  CespereCure  reel-ges  effects  end  nonequilibrium 
thersodynesics  ceused  by  Che  high  teapereCures,  end  perheps  low 
densities  which  cen  be  encountered  in  precCicel  flow  fields. 
Recent  work  in  this  erees  hes  seen  the  incorporeCion  of  resi¬ 
gns  equilibrium  Cheraod3maaics  in  problems  with  two  independent 
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variables  in  an  alaiost  routine  isanner.  Initial  work  has  begun 

on  including  nonequilibrluai  thermodynamics  in  tvc  independent 

31  -}2 

variable  problems.  The  papers  by  Chu  and  Broer  discuss  the 

'  formulation  of  the  problem,  idiile  the  papers  by  Sedney  et.al.^^*^^, 

35  36 

^  Capiaux  and  Washington  ,  and  Wood  et.al.  also  give  some  initial 
results  of  numerical  calculations.  To  the  author's  knowledge,  no 
work  has  been  published  on  problems  involving  more  than  two  in¬ 
dependent  variables  together  with  nonequilibrium  thermodynamics. 

The  object  of  thiv  work,  therefore,  is  to  formulate  the 
method  of  characteristics  for  general  three-dimensional  unsteady 
flow  problems,  allowing  for  the  possibility  of  nonequilibrium 
thermodynamics.  It  is  also  the  purpose  of  this  work  to  write  and 
check  a  working  computer  program  which  applies  the  method  to 
practical  flow  problems.  It  may  be  impossible  or  iaq>ractical  at 
this  time  to  write  a  completely  general  program,  but  only  those 
compromises  and  approximations  which  are  absolutely  necessary  to 
the  writing  of  a  working  program  will  be  stade.  In  writing  the 
program,  it  will  always  be  borne  in  mind  that,  even  though  some 
compromises  must  be  made,  they  will  eventually  be  removed  at  a 
future  time.  In  this  way,  the  compromises  and  approximations 
can  be  introduced  in  such  a  way  that  they  can  be  more  easily 
removed  later. 
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CHAPTER  2 


THE  EQUATIONS  FOR  THE  GENERAL  METHOD 

The  quesl-llncer ,  first-order  partlel  differential  equa¬ 
tions  for  the  flow  of  an  inviscid  nonequilibriua  gas  are  set 
forth  in  this  chapter.  A  brief  review  of  the  theory  of 
characteristics  is  presented  in  Appendix  A.  The  theory  is 
used  in  this  chapter  to  derive  the  characteristic  and  cow- 
patibility  equations  which  are  required  for  the  method  of 
characteristics • 

2.1  The  Equations  of  Change 

The  aquations  which  govern  the  flow  of  a  nonequilibrium 
gas  are  the  conservation  laws,  together  with  the  equation  of 
state.  In  cartesian  coordinates  (x,  y,  s,  t)  they  have  the 
following  fora. 

continuity 

momentum  conservation  in  x-direction 

k-t  ix  ill  di  r  dx  (2.1b) 


(2.U) 
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■OBcntua  conservation  in  y-dlrectlon 


ir  .  /.r  ^  V-  i 

^!t  ilK  diS~^ 


(2.1c) 


■oaentua  conservation  in  z-direction 


^  M  V'^\^4  i-  ifiL 

d)C  ^ 


(2. Id) 


energy  conservation 

6fi .  u  ^  ^  4-  l^^dh 

it  dK  ^3 


E  +  vi£ 

71  it  5 

X  4a 

dBj 

continuity  of  nth  species 

W  ^  V"  ^r\  «  i  (7^  0 

dt  6k  ^3  di  y 

(n  *  2^ . . . .  j 

state  equation 

^  ^  ••  •  •  ?  Ck  ) 


(2.1e) 


(2. If) 


(2.1g) 


Where  in  general 

0'h’‘O’r,(t>,fjCi, - ,C«)  (n.l,2,...jN)  (2-ih) 

is  a  given  function  detenined  froa  the  cheaiatry  of  the  gas 
being  considered.  The  chemical  source  function^  ^n*  *  com¬ 

plicated  function  of  the  chemical  reactions  and  reaction  rates 

39 

being  considered  These  are  2N  +  6  equations  in  2N  *•*  6 
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unknowni  and  thus  form  a  coaplate  act. 

It  should  be  noted  that  the  species  continuity  equation 
(2 .If)  ehich  is  soMetlnes  called  the  rate  equation  is  really 
■ore  general  than  it  night  appear  above.  Vibrational  relasa- 
tlon  and  eleaent  ionisation  are  also  described  by  equations 
vith  the  saae  fon  at  (2. If).  These  equations  are  obtained  by 
utilising  vibrational  energy  or  electron  concentration  instead 
of  the  species  mass  fractions  and  by  using  the  proper  source 
function,  (S  . 

In  order  to  eliainate  the  two  algebraic  equations  (2.1g-h) 
they  auit  be  substituted  into  equations  (2.1e)  and  (2. If)  re¬ 
spectively.  (2.1e)  then  takes  the  form 


■here 


Li  i.  +  uL  +  vL  +  w'L- 
Dt  it  dK  da 

and  it  is  understood  that  all  other  theraodynaaic  variables 
are  held  constant  when  the  enthalpy,  h,  is  partially  dif¬ 
ferentiated  with  respect  to  any  one  theraodynaaic  variable,  a 
is  the  speed  of  sound  in  the  gas  with  the  aass  fractions  of  the 
species  held  constant,  or  in  other  words,  the  frosen  speed  of 
sound.  It  is  given  by 


-[ihl//fihl  - 1] 


(2.1J) 
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This  Is  ch«  cosiplcts  set  of  partial  differential  that  is  to  be 
solved. 

2.2  The  Characteristic  Equation 

As  indicated  in  Appendix  A,  the  characteristic  equation  is 
found  by  first  introducing  a  coordinate  transforaation  to 
^2*  ^3*  ^4^  coordinates.  It  is  also  useful  to  introduce  a  new 
independent  variable  vith  the  fora 

t'  •  Ut  (2.2) 

where  U  is  an  arbitrary  constant  velocity  used  to  give  t'  the 
diaension  of  a  space  coordinate.  This  allows  the  tlae  coordinate 
and  the  space  coordinates  to  be  treated  siailarly.  Choose  the 
new  coordinates  such  that 

-  constant  (2.3) 

is  a  characteristic  hyper surf ace*.  The  characteristic  equa¬ 

tions  for  the  equations  of  change  (2.1)  is  obtained  by  evalua¬ 
ting  the  deterainant  given  in  equation  (A. 5),  Appendix  A.  The 
result  is 


*  The  hypersurface  in  this  case  is  a  three -diaensional 
aanifold  in  a  four-diaensional  space. 
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Froa  (2.4)  it  can  b«  •••&  that  th«r«  two  Mts  of  roal 
oharactaritcic  hyparaurfaoaa  corratpondlng  to 


n/  1/  ^ 
6^  di  dt' 


O 


(2.5) 


Equation  (2.5)  la  tha  dot  product  of  tha  vactor  ^  (a 
ganarallaad  valoelty  vactor) 


4-ur^C+Ui  (2.7) 

ulMra  "S,  Ic,  t  ara  tha  unit  vactorf  along  tha  z,  y,  m,  t* 
ooordlnata  azat  raipactlvaly ,  with  tha  vactor 

7A<  =  MiT  T+  +  ^,  7  <2.«) 

P  hi  di  it 

7  0]^  !•  nomal  to  tha  charactaristlc  hyparaurf aca .  ^  la  tangant 

to  tha  partlcla  llna  or  world  Ilna,  at  It  la  callad  in  ralatlvlty, 
(aaa  lafaranea  38,  p.  5,  and  p.  114).  Tharafora,  hyparturfacat 
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made  up  of  particle  lines  are  characteristic.  Note  that  the 
particle  line  is  a  (N  +  3) -fold  characteristic  because  the 
corresponding  term  in  (2.4)  is  raised  to  the  N  3  power.  This 
property  is  used  to  great  advantage  in  the  method  of  character¬ 
istics. 

The  second  set  of  characteristic  hypersurfaces  is  given  by 
(2.6).  This  is  a  quadratic  hypersurface,  and  is  a  grAeralisation 
of  the  Hach  conoid  which  appears  in  three-dimensional  steady 
supersonic  flow  and  two-dimensional  unsteady  flow  (see  Refer¬ 
ence  13,  pp.  642-657)  and  is  generally  referred  to  as  the  Mach 
hyperconoid.  If  the  equation  of  the  Mach  hyperconoid  is  %rrltten 
in  the  form 


^1  *  ^1 

the  differential  equation  which  defines  the  hyperconoid  is  (2.6). 
Of  course,  this  equation  can  not  be  solved  in  general  for 
but  at  a  particular  point  (x^,  y^,  z^,  t'^)  the  hyperconoid  is 
tangent  to  the  local  Mach  hypercone  and  is  just  an  envelope  of 
sound  waves  emitted  from  the  point.  The  equation  of  the  Mach 
hypercone  is  easily  derived  from  the  geometry  of  sound  waves 
assuming  a  uniform  steady  flow  with  conditions  the  same  every¬ 
where  as  at  the  point  being^considered ,  and  is  given  by 

U  -i-S" 

+  V ■  ^0  (2.9) 


17 


The  Mach  hyperconc  can  be  used  as  a  local  approximation  to  the 
Mach  hyparconold. 


2,3  The  Compatibility  Equations 

The  compatibility  equations  are  obtained  by  forming  a 
linear  combination  of  the  liquations  of  change  (2.1)  and  deter > 
mining  the  linear  coefficients  such  that  partial  derivatives  in 
a  direction  normal  to  the  characteristic  surface  are  eliminated. 
The  equations  for  the  linear  coefficients  and  compatibility  con¬ 
dition  are  derived  in  Appendix  A  and  given  by  (A. 9)  and  (A.  10). 
There  are  two  sets  of  compatibility  equations,  one  for  each  type 
of  characteristic  hypersurface.  Actually,  there  is  only  one 
aquation  corresponding  to  the  Mach  hyperconoid,  i«hile  there  are 
II  ‘f  3  equations  for  the  particle  line.  See  Reference  38,  pp. 
106-107  for  a  discussion  of  tho  number  of  compatibility  equa¬ 
tions  mhich  correspond  to  a  particular  characteristic  hyper¬ 
surface. 


The  linear  combination  of  the  equations  of  change  is  given 
by  (A.  7)  and  has  the  form 


(2.10) 


For  the  given  equations  of  change,  (2.1),  the  compatibility 
aquation  for  the  Mach  hyparconold  has  the  folloving  terms: 


r  V?. 


(2.1U) 


IS 


Aijj if" 74^,) 

(2.11b)  . 

(2.11c) 

(2. lid) 

(2. lie) 

(2. Ilf) 

fn  =•  1^ ,  0  hjJ 


where 


end 


elso 


?.7/S.,= 

OX 

u/  .  U 

61  6t' 

(2.12) 

%dien  m  ■  1 

D„ 

when  m  "  2 

%rhen  m  ■  3 

when  m  "  4 

S  ~  4  2 

r  h.t 

(2.13) 

These  equeCions  correspond  to  rquetlons  (A.8)  in  Appendix  A. 
The  linear  coefficients,  ^  ^  ,  are  determined  by  setting 
equal  to  zero  and  solving  the  resulting  set  of  equations. 
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This  •llminates  derivatives  in  a  direction  noraal  to  the 
characteristic  hypersurface.  The  result  is 


=  -f 

-  -f  '^/ fv*  7^) 

K/K  =  -f 

Ar/ 1 


(2.14a) 

(2.14b) 

(2.14c) 

(2.14d) 


>  fh»l,  ...‘iN)  (2.14e) 


ifhen  equations  (2.11)  through  (2.14)  are  substituted  into  (2.10), 
the  folloiring  coapatibility  equation  for  the  Mach  hyperconoid 
results. 


(2.15) 
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where 

"  constant 

is  the  characteristic  Msch  hyperconoid.  Note  that  derivatives 
of  ^  and  0^  do  not  appear  in  the  equation. 

As  mentioned  above,  there  are  N  3  compatibility  equations 
for  the  particle  line,  but  only  N  +  1  of  these  equations  arc 
needed  for  the  method  of  characteristics.  These  could  be  ob¬ 
tained  in  a  manner  similar  to  the  procedure  Just  used  to  obtain 
equation  (2.15).  However,  because  the  equations  sought  do  not 
have  derivatives  in  a  direction  normal  to  the  particle  line  and 
thus  involve  derivatives  along  the  particle  line  only,  it  is 
simpler  to  use  this  property  to  choose  them  by  inspection  from 
the  original  partial  differential  equations. 

The  substantial  derivative 

Ot  dt' 

is  the  derivative  of  a  dependent  variable  following  a  fluid 
particle,  so  if  ^  is  taken  as  a  coordinate  along  the  particle  line 
in  the  (x,  y,  z,  t')  space,  equations  (2.1£)  can  be  written 

(2.16) 

where 

</  =  >/u*’+  +  («/*■  +  U*' 
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Thes«  If  of  the  coapatibility  equatioas  along  the  particle 
line.  One  additional  relation  it  needed.  The  energy  equation, 
which  can  be  written  as  (2. If)  or  as  (2.1i),  meets  this  need. 

The  energy  equation  can  be  written  in  a  somewhat  simpler  form  if 
the  specific  entropy  is  Introduced.  The  entropy  is  a  function 
of  the  thermodynamic  variables  and  the  mass  fractions  of  species 

S  s  S  (f>,  j5,  Cl, . Cm)  (2.17) 


SO 

where  the  subscripts  indicate  the  variables  held  constant  in 
the  partial  differentiations.  Substituting  for  Dp/Dt  in  (2.11) 
and  dtlliaing  another  fora  of  the  expression  for  the  froaen 
speed  of  sound 


the  following  is  obtained 


Sttbbtituting  for  Dc^/Dt  from  (2. If)  and  introducing  the  ^  co¬ 
ordinate  gives 
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BquAtion  (2.l8c)  can  alio  be  written  in  temui  of  the  chemical 
potential  of  the  epeclei,  and  the  temperature,  T,  (ate 

Reference  39#  equation  (21)  ) . 

U 

(2.l8d) 

n.l 

Equations  (2.l6)  and  (2.l8c)  are  the  required  compatibility 
relations  for  the  particle  line. 

2,4  The  Nonequlllbrlum  Terms 

As  was  mentioned  above  In  Section  2.1,  the  chemical  source 
function  Is  a  complicated  function  of  the  exact  chemical 
reactions  being  considered  and  their  reaction  rates,  llie 
reaction  ratesi.  In  turn,  can  be  a  function  of  temperature. 

The  thermodynamic  functions 

h  -  h  (p,  ^  - -  Cj^) 

S  -  8  (p,  ^  ,  C^,  - -  Cjj) 

can  also  be  conpllcated  functions  of  the  mass  fractions. 
Specific  foirnis  of  these  various  functions  are  not  considered 
here  as  they  vary  with  the  type  of  reactions  being  considered. 
The  general  forms  of  the  reaction  rates  In  certain  simple 
cases  can  be  obtained  analytically,  but  In  practice  for  most 
reactions,  they  must  be  determined  experimentally.  Because 
of  this,  many  of  the  functional  relationships  are  unknown,  or 
at  least  uncertain  at  the  present  time.  Therefore,  only 
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g«n«ral  functional  formi  are  uaed  here.  Thlf  in  no  way 
llmlti  the  results  of  the  method  of  characteristics  formu¬ 
lated  here.  When  a  specific  case  is  being  considered,  it 
is  merely  a  matter  of  substituting  the  proper  functional 
relations  in  the  appropriate  equations  in  order  to  obtain 
the  goYeming  partial  differential  equations.  Note  that 
we  are  referring  to  the  partial  differential  equations  here. 
The  solution  of  the  equivalent  finite  difference  equations 
is  by  no  means  slnqple,  as  is  pointed  out  in  Section  3*2 
of  the  next  chapter. 
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CHAPTER  3 


THE  ELEMEWTS  OF  THE  INTEGRATION  PROCEDURE 

In  this  chapter,  the  partial  differential  equations  devel¬ 
oped  in  Chapter  2  are  written  in  finite  difference  foxv  and  the 
numerical  procedures  for  solving  the  resulting  equations  are 
presented.  The  process  involves  choosing  a  finite  difference 
network  and  then  developing  the  finite  difference  aquations 
which  correspond  to  it. 

3.1  The  Choice  of  the  Metwork  Configuration 

The  method  of  characteristics  for  two  independent  variables 
has  been  used  extensively  and  involves  rather  simple  geometrical 
principles.  In  considering  more  than  two  independent  variables, 
the  geometry  of  the  characteristic  hyperturfaces  becosMS  more 
complicated  and  some  fundamental  changes  have  to  be  made  to  ex¬ 
tend  the  method  to  onxlti-variable  problems.  First,  we  shall 
consider  the  extensions  required  to  go  from  two  to  three  in¬ 
dependent  variables  after  which  the  extensions  to  four  variables 
can  be  more  easily  determined. 

Consider  the  geometry  of  the  characteristic  lines  for  two 
independent  variables  as  shown  in  Figure  1.  If  the  initial 
data  are  given  on  the  line  P2^P2  wt  points  P^^  and  P2»  there  is 
only  one  possible  characteristic  net  available  to  determine 
point  P^*  This  is  because  only  two  characteristic  lines  pass 
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through  each  point  and  the  only  possible  adjacent  intersection 
on  the  right  side  on  the  initial  data  line  is  at  point  P^*  In¬ 
troducing  the  streamline  or  particle  line,  which  is  also  a 
characteristic  line,  docs  not  allow  other  network  configuration 
cither.  Hence,  in  obtaining  the  solution  by  marching  to  the 
right  from  the  initial  data  line  P]^P2  only  one  possible 

characteristic  net  which  can  be  used. 

This  is  not  the  case  in  problems  involving  three  independent 
variables.  When  the  extra  variable  is  added,  an  extra  "degree 
of  freedom"  enters  which  allows  a  choice  in  the  net  configuration. 
This  can  best  be  seen  by  referring  to  Figure  2.  One  sheet  of  the 
Mach  conoid  (e.g.  in  the  downstream  or  positive  time  direction) 
through  a  point  P  is  depicted.  This  conoid  is  a  special  case  of 
the  second  characteristic  surface  given  by  (2.6),  but  with  one 
of  the  independent  variables  omitted.  This  Mach  conoid  is  Just 
the  envelope  of  all  the  local  infinitesimal  Mach  cones.  As 
shown  in  Figure  3,  the  local  Mach  cones  are  tangent  to  the  Mach 
conoid  at  every  point  on  its  surface.  A  line  on  the  Mach  conoid 
which  is  everywhere  tangent  to  the  generators  of  the  local  Mach 
cones,  is  termed  bicharacteristic  and  is  shown  in  Figure  2.  It 
is  seen  that  there  is  a  single  parameter  family  of  bicharacter¬ 
istics  on  the  conoid  surface.  The  parameter  defining  a  bi- 

eharacteristic  might  be  chosen  as  the  angle  9  shown  in  Figure 

I 

2. 

Hext,  it  is  useful  to  introduce  the  concept  of  character¬ 
istic  surfaces  corresponding  to  a  curve.  Consider  a  specified 
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curve  PQ  8ho%in  in  Figure  4  together  with  the  conoids  correspond¬ 
ing  to  each  point  on  the  curve.  The  envelopes  of  the  upper  and 
lower  portions  of  the  conoids  are  the  characteristic  surfaces 
corresponding  to  the  curve  PQ.  Note  that  the  lines  of  tangency 
of  these  surfaces  with  the  conoids  are  the  bicharacteristics. 

Now,  returning  to  Figure  2,  consider  the  normal  vector 
N  which  is  normal  to  the  conoid  surface  at  P  and  hence,  is 
noLmal  to  the  tangent  to  the  bicharacteristic  at  P.  is  also 
normal  to  one  of  two  characteristic  surfaces  corresponding  to  a 
line  through  P  and  normal  to  N.  Thus  for  each  normal  vector  and 
bicharacteristic,  there  is  at  least  one  characteristic  surface 
corresponding  to  a  line  through  P.  There  is  at  least  a  single 
parameter  family  of  such  surfaces  through  every  point,  P. 

Finite  difference  net  configurations  can  be  made  up  of  these 
surfaces  and  because  of  the  availability  of  an  infinite  number 
of  such  surfaces  through  each  point,  there  is  an  additional 
"degree  of  freedom"  in  the  choice  of  the  net. 

It  is  a  simple  extension  to  reason  that  when  a  fourth 
independent  variable  is  included,  a  second  degree  of  freedom  in 
the  choice  of  the  net  configuration  is  added.  This  is  because, 
corresponding  to  the  characteristic  surfaces,  characteristic 
hypersurfaces  (three-dimensional  manifolds  in  the  four-space) 
could  be  obtained  at  each  point,  and  a  two  parameter  family  of 
such  hypersurfaces  would  be  available.  There  would  then  be  a 
double  infinity  of  hyper surf aces  to  choose  from  in  making  up  a 
net  configuration. 
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With  this  understanding  of  the  generalisation  of  the  finite 

difference  net  to  probleais  involving  three  end  four  independent 

veriables,  the  choice  of  a  specific  net  configuretion  can  now  be 

attesipted.  Again,  consider  the  three  independent  variable  case 

before  trying  the  nore  coapliceted  four  variable  problem.  Various 
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network  configurations  can  be  proposed.  Fowell  has  reviewed 
several  networks  which  have  been  proposed  and  in  some  cases, 
utilised  by  various  people.  We  will  not  go  into  the  details 
of  the  networks  here,  but  will  give  a  very  brief  description  of 
each  and  compare  their  advantages  and  disadvantages.  For  further 
details,  the  reader  is  referred  to  Fowell  and  the  original  papers 
which  first  described  each  network. 

Thomhill^^  first  proposed  two  network  configurations  which 
Fowell  has  termed  the  tetrahedral  characteristic  line  network  and 
the  tetrahedral  characteristic  surface  network.  The  nomenclature 
introduced  by  Fowell  is  adopted  here.  These  networks  are  shown 
in  Figures  5  and  7.  A  minimum  of  three  points  is  required  in 
the  initial  surface  to  calculate  conditions  at  a  new  point  off 
the  surfece.  In  the  tetrahedral  characteristic  line  network, 
three  Mach  cones  through  the  three  initial  points,  and 

P^,  are  utilised  to  form  the  net.  A  mutual  intersection  point 
of  the  three  cones  is  the  new  point  at  which  the  flow  proper¬ 
ties  are  to  be  calculated.  The  generetors  of  the  cones  which 
pass  through  the  intersection  point  are  approximations  to  bi- 
characteristics.  The  competibility  equations  can  be  written  in 
finite  difference  form  along  these  bicharacteristir.s  and  then 
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utlllEcd  to  dttoninc  the  flow  properties  et  the  new  point. 

If  this  network  is  considered  frow  the  standpoint  of  nuwericel 
stability,  however,  it  is  found  to  be  unstable  (see  Appendix 
B).  For  this  reason,  this  writer  proposes  the  modified 
tetrahedral  characteristic  line  network  shown  in  Figure  6. 

In  the  Bodified  network,  it  is  proposed  that  the  triangle 
foraed  by  the  three  initial  points  be  inscribed  with  a  circle, 
and  the  points  of  tangency  of  the  circle  with  the  triangle,  F^2» 
P23»  end  be  utilized  in  exactly  the  saae  way  as  the  initiAl 

points,  P^,  P2*  end  P^,  for  the  unmodified  network.  The  prop¬ 
erties  et  the  points  ^^2*  ^23*  ^31  ^  obtained  by 

linear  interpolation  between  the  initial  points  or  a  higher  order 
interpolation  scheae  might  be  utilized  in  the  initial  surface. 

The  addition  of  the  procedure  of  inscribing  the  circle  within 
the  triangle  assures  numerical  stability  for  this  network  (see 
Appendix  B). 

✓  ^ 

The  tetrahedral  characteristic  surface  network  also  uses 
the  ■Inimum  of  three  points  in.  the  initial  surface  to  form  the 
net.  The  lines  passing  through  each  pair  of  points  are  used  to 
define  characteristic  surfaces.  The  mutual  intersection  point 
of  the  three  inward  leaning  characteristic  surfaces  is  utilized 
as  the  new  point,  P^.  The  characteristic  conoid  through  P^ 
opening  toward  the  initial  surface  is  tangent  to  the  character¬ 
istic  surfaces  along  bicharacteristics  as  shown  in  Figure  7. 

The  points  F^2*  ^23*  ^31  base  of  the  bicharacteristics 

in  the  initial  surface  are  utilized  in  the  finite  difference 
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•quaCioDS  and  the  flow  properties  at  these  points  must  be  deter- 
■Ined  by  interpolation  in  the  initial  surface. 

26  18 

Ferrari  and  Sauer  have  proposed  a  network  made  up  of 
characteristic  surfaces  and  two  orthogonal  families  of  reference 
planes  as  shown  in  Figure  8.  This  network  of  intersections  of 
reference  planes  with  characteristic  surfaces  utilises  two 
families  of  coordinate  planes  (e.g.  x  *  constant  and  s  -  constant) 
as  depicted  in  Figure  8.  The  points  P^,  ?2f  F3  ond  in  the 
initial  surface  are  used  to  calculate  the  flow  at  the  new  points 
P^  and  P^  which  are  located  at  the  intersections  of  a  "left 
running"  characteristic  surface  from  P]^F2*  *  "right  running" 
characteristic  surface  from  F3F4  ond  a  specified  x  ■  K  co¬ 
ordinate  plane.  The  flow  properties  at  P^,  for  instance,  are 
obtained  from  equations  for  the  variations  of  flow  properties 
along  the  lines  P^P^  and  P3F5  in  the  characteristic  surface. 

The  third  equation  is  obtained  by  considering  the  variation  in 
the  initial  surface  (along  P|^P2  or  P3F4)  nnd  relating  it  to  the 
variations  along  P]^F3  ond  ^3^3*  The  flow  properties  at  P^  are 
obtained  similarly,  but  note  that  the  base  point  P^  will,  in 
general,  have  to  be  moved  in  an  iteration  to  Pg  in  order  to 
cause  Pg  to  fall  on  the  x  *  R  coordinate  plane.  This  will  re¬ 
quire  interpolation  or  extrapolation  in  the  initial  surface  in 
order  to  find  the  flow  properties  at  Pg. 

The  prismatic  network  of  characteristic  surfaces  was  intro¬ 
duced  by  Cobum  and  Dolph^^  and  refined  by  Holt^^.  It  is  shown 
in  Figure  9.  This  network  is  similar  to  the  network  of 
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intersections  of  reference  planes  with  characteristic  surfaces 
except  that  bicharacteristics  are  used  in  integrating  in  the 
characteristic  surface  rather  than  the  intersections  «rith 
orthogonal  coordinate  planes.  This  network  also  requires  that 
the  flow  be  known  on  a  surface  other  than  the  initial  surface 
such  as  surface  P2^4^6  figure.  This  might  be  a  plane 

of  symmetry,  but  it  seems  to  limit  the  useful  application  of 
the  network.  The  flow  is  determined  at  the  new  point  by 
solving  equations  «n:itten  along  the  lines  ^6^5* 

28 

Butler  has  proposed  a  network  and  method  which  is  not 
precisely  the  method  of  characteristics,  but  is  more  of  a  com¬ 
bination  of  the  standard  finite  difference  approach  and  the 
method  of  characteristics  in  the  spirit  of  the  method  of 
Hartree^  mentioned  in  Chapter  1.  Foi#ell  termed  this  the 
pentahedral  bicharacteristic  line  network.  As  shown  in  Figure 
10,  the  coordinates  of  the  new  point  Pq  are  chosen, and  the 
Kach  cone  is  projected  back  toward  the  initial  surface.  Four 
bicharacteristics,  one  more  than  the  minimum  number  required, 
are  used  in  determining  the  flow  properties.  The  addition  of 
the  fourth  equation  allows  the  elimination  of  the  partial  deriva 
tives  of  the  flow  properties  at  the  new  point.  The  flow  pro¬ 
perties  at  base  points,  P^,  ?29  P3  and  P^,  which  lie  in  the 
initial  surface,  must  be  obtained  by  interpolation.  Note  that 
Che  location  of  the  bicharacteristics  is  not  arbitrary,  but  is 
determined  by  the  condition  for  the  elimination  of  the  partial 
derivatives  so  that  the  simplification  of  eliminating  some 
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dtrivativct  it  traded  for  the  complication  in  determining  the 
location  of  the  bicharacteri sties  and  the  addition  of  an  extra 
aquation. 

24  40 

Sauar  and  Schaetz  have  proposed  methods  which  utilize 
what  they  call  "near  Characteristics",  but  %rhich  are  very 
similar  bo  the  prismatic  network  of  characteristic  surfaces  and 
tha  natwork  of  intersections  of  reference  planes  with  character¬ 
istic  surfaces  discussed  above.  It  is  felt  that  the  near 
characteristic  networks  have  advantages  and  disadvantages 
similar  to  those  discussed  above  for  the  prissmtic  and  reference 
plana  networks. 

Fowall  has  discussed  in  detail,  the  advantages  and  draw¬ 
backs  of  each  of  the  networks  outlined  above.  We  only  summarize 
hare,  his  discussion  adding  a  point  about  stability  of  which  he 
was  not  aware.  In  general,  three  points  can  be  used  to  compare 
tha  networks.  They  are:  (1)  Is  interpolation  required  in  the 
initial  surface?  This  is  to  be  avoided  because  it  can  introduce 
inaccuracy.  (2)  Do  the  base  points  in  the  initial  surface  move 
about  as  tha  solution  is  iterated  to  higher  order  in  the  step 
size?  If  the  points  do  move,  interpolation  is  required.  (3)  Are 
the  equations  to  be  solved,  written  along  bicharacteristic  curves? 
The  formal  proofs  of  existence  and  uniqueness  given  by  Titt^^ 
apply  only  for  integration  along  bicharacteristics. 

Interpolation  is  required  for  all  the  networks  except  the 
tetrahedral  characteristic  line  network.  Unfortunately,  it  is 
numerically  unstable.  The  modified  tetrahedral  characteristic 


line  network  requires  only  one  interpolation  so  that  it  is  the 
next  best.  The  others  require  considerably  more  interpolation. 

The  base  points  must  be  moved  for  the  last  three  networks 
outlined  above  and  sho%m  in  Figures  8,  9,  and  10.  They  are 
fixed  for  the  first  three  networks. 

Finally,  on  the  third  point,  bicharacteristics  are  used 

in  all  but  the  network  of  intersections  of  reference  planes 

with  characteristic  surfaces.  It  should  be  noted,  however, 
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that  Morretti  et.al.  have  used  this  net  and  have  been  success* 
ful  in  obtaining  practical  results  so  that  the  argument  for 
using  bicharacteristics  is  not  clearly  established.  It  should 
also  be  noted  that  the  new  points  are  obtained  on  equally  spaced 
coordinate  surfaces  using  this  net  so  that  boundary  conditions 
can  he  more  easily  handled  in  some  problems. 

There  are  other  advantages  and  disadvantages  of  these  net* 
work  configurations  that  cannot  be  presented  here  because  of 
space  limitations,  but  the  major  points  have  been  discussed 
above.  Fotrell  felt  that  the  tetrahedral  characteristic  line 
network  was  the  most  favorable  and,  at  least  in  a  relative  sense, 
that  its  advantages  outweighed  its  disadvantages.  With  the  dis* 
covery  of  the  numerical  instability  in  this  network,  this  writer 
believes  that  the  modified  tetrahedral  characteristic  line  net* 
work  is  now  the  favored  approach.  Admittedly,  this  sort  of 
reasoning  is  very  qualitative,  but  it  appears  to  be  the  best  that 
can  be  done  short  of  programming  all  of  the  networks  and  com* 
paring  the  results  they  give  for  a  test  problem.  Of  course,  the 
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latter  suggestion  is  out  of  the  question  practically. 

Now  that  we  have  seen  the  difficulty  in  choosing  a  net¬ 
work  configuration  for  problems  involving  three  independent 
variables,  it  should  be  apparent  that  choosing  a  network  for 
a  four  variable  problem  is  even  more  difficult.  First  of  all, 
our  inability  to  "picture"  the  geometry  of  a  four  variable 
space  coiq>licate8  the  choice.  This  does  not  make  the  task 
Impossible,  however,  because  the  process  can  still  be  approached 
on  a  mathematical  basis.  It  might  be  expected  that  a  general¬ 
isation  of  the  modified  tetrahedral  characteristic  line  network 
would  also  be  a  logical,  advantageous  network  to  use  in  four 
independent  variables.  Proceeding  with  this  point  of  view,  we 
can  devise  the  generalised  network. 

Consider  four  points,  P^,  ?2»  P3  luid  P^  lying  in  the 
initial  hyper  surf  ace.  Four  is  the  minimum  number  of  points  re¬ 
quired  to  determine  the  flow  properties  at  the  new  point  P^. 
First,  corresponding  to  the  circle  inscribed  within  the 
triangle  in  the  three-variable  networks,  we  inscribe  a  sphere 
within  the  tetrahedron  formed  by  the  initial  points  in  the 
initial  hypersurf ace .  The  four  points  of  tangency  of  the  sphere 

with  the  tetrahedron  labeled  ^2^23*  ^124*  ^134  ^234 

"representation"  in  Figure  11(a)  are  used  as  the  base  points 
for  four  Mach  hypercanes  shown  in  Figure  11(b).  The  flow 
properties  at  these  four  points  can  be  obtained  by  interpola¬ 
tion  in  the  initial  hypersurface.  The  equation  of  the  Mach 
hypercones  (2.9)  is  a  second  order  algebraic  equation  in  four 
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unknoims.  This  equation  written  for  each  of  the  four  inter* 
mediate  base  points  (i  ■  1,  2,  3,  4)  gives  four  equations  in 
four  unknowns  which  can  be  solved  for  the  coordinates  of  the 
Intersection  point,  P^.  The  compatibility  equation  (2.15) 
can  be  written  in  finite  difference  form  along  the  four  lines 
from  the  intermediate  base  points  to  the  new  intersection 
point.  This  gives  four  equations  to  determine  four  flow  prop¬ 
erties  at  the  new  point,  for  example,  u,  v,  w,  and  p.  The  re¬ 
maining  flow  properties  at  the  new  point,  such  as  the  entropy, 
s,  and  the  species  mass  fractions,  c^,  are  determined  from  the 
compatibility  equations  (2.16)  and  (2.18c)  which  apply  along 
the  particle  line  through  the  new  point  projected  back  to  its 
intersection  «rith  the  initial  hyper surf ace.  The  flow  proper¬ 
ties  at  the  base  of  the  particle  line  in  the  initial  hyper¬ 
surface  must  also  be  obtained  by  interpolation. 

This  completes  the  discussion  of  the  network  configuration 
to  be  used,  but  before  proceeding  to  the  details  of  the  finite 
difference  equations,  it  should  be  pointed  out  that  new  points 
on  various  types  of  boundaries  must  also  be  calculated.  Points 
on  shock  waves,  on  body  surfaces,  on  contact  surfaces,  and  on 
other  discontinuities  which  might  arise  in  flow  fields  of 
interest  must  also  be  calculated.  The  networks  described  above 
are  for  what  might  be  called  field  points,  which  are  in  the 
interior  of  flow  fields  away  from  the  boundaries  and  interior 
discontinuities.  For  the  flow  fields  of  interest  here,  at 
least  two  other  types  of  points  must  be  calculated.  They  are 
points  which  lie  on  a  body  surface  and  points  %diich  lie  on  an 
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•xtsrlor  shock  wmvc  surface.  By  exterior  shock  wave  surface, 
ve  aean  a  surface  beyond  which  the  flow  is  known  or  un¬ 
disturbed.  Interior  sccpndary  shock  waves  will  not  be  con¬ 
sidered  in  this  initial  study,  although  they  can  be  calculated 
with  SOM  added  effort.  The  three  types  of  points  to  be  con¬ 
sidered  here,  are  shown  in  the  three -diaensional  representation 
of  the  four -space  in  Figure  12  fof  the  flow  in  the  nose  region 
of  a  blunt  body  saving  at  a  supersonic  velocity. 

Note  that  in  several  of  the  figures,  three-dinensional 
"representations"  are  used  for  the  four -dimensional  space. 

Thus,  a  two-dimensional  surface  may  be  draim  in  a  figure,  but 
it  may  represent  a  hypersurface  and  will  be  labeled  accordingly 
as  in  Figure  12  where  the  body  and  shock  hypersurface  are  so 
labeled.  Surfaces  depicted  in  the  figures  may  indicate  a  sur¬ 
face  or  they  may  represent  a  hypersurface,  but  lines  depicted 
in  the  figures  always  represent  lines.  In  Figure  12,  the 
initial  hypersurface  is  represented  by  the  hyperplane  t*  -  0 
(the  initial  steady  flow). 

3.2  The  Field  Point  Procedure 

In  this  section,  the  details  of  determining  the  coordinates 
and  flow  properties  at  a  new  point  in  the  field  are  presented. 
The  solution  of  a  new  field  point  requires  the  completion  of 
four  steps.  First,  the  intermediate  base  points  must  be  located 
and  the  flow  properties  at  each  point  detemined.  The  location 
of  the  new  field  point  is  then  detemined  and  as  the  third  step, 
the  flow  properties  at  the  new  point  must  be  obtained.  Finally, 


a  coavcrgant  itaratioa  schoM  ausC  be  establlahad  which  will 
datenilnc  the  location  and  flow  properties  at  the  new  point 
to  second-order  accuracy  in  the  step  size.  The  details  of 
each  of  the  steps  follows. 

Step  1 

The  four  base  points  determine  a  tetrahedron  within  which 
a  sphere  is  to  be  inscribed.  The  tangency  points  are  utilized 
as  intermediate  base  points  from  which  the  Hach  hypercones  are 
projected.  The  flow  properties  at  the  tangency  points  must  be 
obtained  by  interpolation. 

First,  the  calculation  is  simplified  by  transforming  to  a 
set  of  cartesian  coordinates  (x*,  y*,  z*,  t*)  in  which  the  t* 
coordinates  of  the  four  points  are  equal.  This  is  possible 
because  the  four  points  cannot  lie  in  more  than  a  three-dimen¬ 
sional  subspace  of  the  four-space.  In  this  subspace,  the  planes 
which  are  the  faces  of  the  tetrahedron  can  be  written  in  the  form 

Ai  X*  i  Cc  i*  =1  (3.U) 

The  sphere  has  the  equation 

=  <3.ib> 

Finally,  the  nonsal  vector  to  the  sphere  ssist  be  parallel  to 
the  normal  to  each  of  Che  faces  of  the  tetrahedron  at  each  point 
of  tangency.  This  condition  is  given  by  the  following  equations. 
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Bc/|fVfl  =  (y*-ae)/^lNsfhl  U-lj...,^) 


(3.1c) 


(3.W) 


=  (Ar  +  B'  + 


(3.U) 


(3.1f) 


Tht  third  coi^>onent8  of  the  norael  vectors  ere  eutoaetlcelly 
equel  beeeute  of  the  norMllsetlon.  Equetlone  (3.1e-d)  epply 
dt  eech  of  the  points  of  tengency,  end  thus  these  four  eque- 
tlons  written  for  eech  of  the  four  points  of  tengency  give  16 
equetloQS  for  the  sixteen  unknowns,  y*^f  (1  •  1,  2, 

3,  4),  X*Q,  y*Q,  x*Q,  end  1. 

This  set  of  equetlons  cen  be  solved  nuaerlcelly  by  using 

e  Hcwton  iteretloo  schesM  (soaetlaes  teraed  the  Newton-Rephson 
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elgorltfasi,  see  Lepldus  ,  p.  288).  This  schesM  requires  e  first 
estlaete  of  the  roots  of  the  equations.  These  cen  be  taken  as 
the  center  of  eech  face  of  the  tetrahedron  for  the  tengency 
points  ^1*  ^1*  center  of  the  tetrahedron  for  the 

center  of  the  sphere  x*q,  g*Q,  g*Q  and  the  average  distance 
froBi  the  center  of  each  face  to  the  center  of  the  tetrahedron 
for  the  radios  of  the  sphere  R. 

After  the  coordinates  of  the  tengency  points  are  deter¬ 
mined  in  Che  (n*,  y*,  a*,  c*)  coordinates,  the  inverse 


3t 


1 

j 

< 

CransforMtion  is  used  to  return  to  the  (x,  y,  z,  t')  systea. 
The  flow  properties  at  the  tangency  points  must  be  obtained 
by  Interpolation.  The  simplest  procedure  would  be  to  linearly 

'  f 

Interpolate  between  the  three  comer  points  on  each  face  of 
the  tetrahedron.  This  could  be  done  In  the  following  manner. 

Starting  In  the  (x*,  y*,  z*,  t*)  coordinates  for  each 
plane  face,  transform  to  ex'**,  z*^)  coordinates.  Drop  the 

t*  coordinate  and  work  In  the  subspace  because  the  t*  co¬ 
ordinates  of  all  poln::8  being  considered  are  equal.  Consider 
the  face  Figure  11(a).  In  the  transformation  to  the 

(x  ,  y  ,  z  )  system  the  origin  Is  located  at  and  the  x 
axis  aligned  with  the  line  following  equations: 

»[Cx*-X,‘')ces6  -h  S.n^Ccs  (3,2.) 

3  ^in  4) 

-fft*- l^Jzindjim^Jen  f  (3.2b) 

*iC -  fxf-XiJ  sinfijsin  f 

i*  -('/*  ^  (3.2c) 

~ XiJceiG  -I-  (i*-  i*J  sin  iift  ^ 

where 

+«n  f  »  ■  ^j)ics6  -  ♦ 

-  +  (a*-2t)5in6]s:n  <pj 

tan  4>  =■  (.yt  ~^t)/i(X*-)(t)c<>s&-h(i*-2*)'!.,n  6'J 
+.n  e  "  (it-itjytxt-Xi) 


(3.2d) 

(3.2e) 

(3.2f) 
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Intarpolating  in  this  plana  with  f  representing  a  typical  flow 
property,  the  property  at  the  tangency  point  is  given  by 


itl 


^2 


'•15 


L 


ylti 


(3.3) 


The  coordinates  and  flow  properties  for  the  Interaediate 
base  points  have  been  determined  and  thus  the  first  step  of  the 
field  point  solution  is  complete. 


S£se_2 


The  coordinates  of  the  new  field  point  are  obtained  by 
simultaneonsly  solving  the  equations  of  four  Mach  hypercones 
(2.9)  which  eminate  from  the  four  intermediate  base  points 
determined  in  Step  1  above.  The  equations  can  be  written  in 
the  form 


(X-  X t  (a  *  b  f  *  -  *i  )*+  (t'-t;- 1  ’  * 


(3.4«) 


where  (i  *  1,  2,  3,  4),  which  designates  each  of  the  inter¬ 
mediate  base  points.  Tliese  are  four  equations  of  second  degree 
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in  four  unknoims  and,  in  general,  have  16  roots.  Only  the 
closest  root  to  the  initial  hypersurface  in  the  positive  tiae 
direction  is  desired. 


To  solve  (3.4a),  consider  it  rewritten  and  expanded  in 

\(0)  V 


the  following  fora: 

+^(K-0(y- 

•  a  .  .  .  S 

where 

F*'  =  F.'  f’.  t'-\  t“<) 

i/®’ 


(3.4b) 


and  the  superscript  in  parentheses  indicates  an  estiaate  of 
the  root  or  a  step  in  an  iteration  process.  Now,  an  iteration 
aight  be  carried  out  using  the  following  t3rpe  of  equations 


M  .("-I) 

<  =  l(  +  4X 

2  '=  £ 


(3.4c) 


where  for  exaaple 

=.  (k- 

in  (3.4b).  Ax  can  then  be  obtained  by  truncating  (3.4b)  at 
soac  point  and  solving  for  Ax.  A  first  estiaate  of  the  root 
(x^°\  t'^®^)  must  be  available  to  start  the  process. 
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Foifsll  found  that,  for  raasonabla  first  attlnatas,  tha 
•ocond  order  derlvatlvai  in  (3«4b)  nuat  be  included  to  aeeure 
convergence.  Fowell  refers  to  this  as  a  second-order  Nevton 
procedure  eliere  second-order  probably  refers  to  the  inclusion 
of  second-order  derivatives.  Actually,  the  Neirton  algor ithn 
includes  only  the  first  order  derivatives  and  it  accurate  to 
second  order  in  the  error  tern  Ax.  The  inclusion  of  the 
second-order  derivatives  is  similar  to  the  Richmond  iteration 
procedure  (see  Lapidus  ,  p.  292)  and  is  good  to  third  order 
in  Ax. 


When  the  second-order  derivatives  are  included,  Fowell 
suggests  iterating  the  following  equations  to  obtain  the 
correction  terms. 


(3.4d) 


o 


Values  of  sero  can  be  taken  as  first  estimates  of  the  cor¬ 
rections.  These  are  substituted  into  the  terms  in  brackets 
and  the  resulting  linear  equations  can  be  solved  for  a  new 
value  of  the  corrections.  These  corrections  are  then  sub¬ 
stituted  into  the  bracketed  terns  and  the  process  repeated 
until  convergence  is  obtained.  The  corrections  are  then 
used  in  (3.4c)  to  obtain  new  values  of  the  coordinates. 
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These  new  coordlnetes  ere  in  turn  used  to  re-eveluete  end 
its  derivetlvcs  in  (3.4d),  end  the  entire  process  is  repeeted. 
Thus,  the  procedure  involves  e  double  iteretion. 


This  procedure  suggested  by  Fowell  wes  utilized  in  the 
celculetions  cerried  out  for  this  study,  but  it  should  be 
pointed  out  thet,  by  following  Richmond's  iteretion  process 
more  closely,  one  iteretion  required  in  the  ebove  double  itere¬ 
tion  cen  be  elimineted.  If  (3.4b)  is  trunceted  efter  the  first 
order  derivetives,  we  obtein 


(3.4e) 


These  equetions  cen  be  solved  for  the  corrections,  end  used  to 
eliminete  the  corrections  eppeering  in  the  breckets  in  (3.4d). 

The  corrections  outside  the  breckets  ere  reteined.  These  new 
equetions  cen  then  be  explicitly  solved  for  the  corrections  so 
thet  iteretion  is  unnecessery.  This  second  procedure  is  e 
little  less  eccurete  in  determining  the  correction  but  eliainetes 
one  iteretion,  though  et  the  expense  of  introducing  more  com- 
pliceted  equetions.  This  mey  be  e  more  efficient  procedure, 
but  it  is  very  difficult  to  determine  definitely  short  of  pro- 
grenming  end  compering  the  two  schemes. 

Following  Fowell,  we  obtein  e  first  estimete  of  the 

/* 

desired  root  of  <3.4e)  with  the  following  reesoning.  Consider 
ell  four  hypercones  to  be  identicel  end  defined  by  properties 
which  ere  evereges  of  the  properties  of  the  ectuel  hypercones. 
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Froa  a  point  an  equal  distance  from  the  four  intermediate  base 
points  lying  in  the  initial  hyperplane,  project  a  line  in  the 
direction  of  the  average  generalized  velocity, 

The  first  estimate  of  the  location  of  the  new  field  point  lies 
out  along  this  line  in  the  positive  time  direction  at  a  distance 
equal  to  d  cotyt/  where  yU  is  an  average  Mach  angle  for  the  cones 
and  d  is  the  distance  from  one  intermediate  base  point  to  the 
point  located  an  equal  distance  from  the  four  intermediate  base 
points.  This  reasoning  results  in  the  following  equations  for 
the  first  estimates  of  the  coordinates  of  the  new  field  point. 


d 

2^’=  ii&'-d 

aave 

where,  for  example^ 


-#• 


y- 


(3.4f) 


4-  he 
^  t'c 


(«i  +  Wi  ^ 


and  the  numbered  subscripts  indicate  each  of  the  intermediate 
base  points.  Also 

d  =  -ti/J* 

the  c  subscripts  indicate  the  coordinates  of  the  point  lying 
in  the  initial  hyperplane  an  equal  distance  from  the  inter¬ 
mediate  base  points.  These  are  determined  from  the  following 
equations 
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(it-Xi)V(yc-ai)'  +  f2t'  2i)%  (t'c  - 

-GCt-xO*- 

~  (Xe-O*"  +(34-  aif+fZt-ijArtfi-  t s)*" 

+  Cijc-yi)%-  ( <de  ~  Iji)^  (-t'c -t'l)^ 

AiCe  +  Bye  +Cic  +  Dtj  s  1 
where  the  lest  equetlon  is  the  equetion  of  the  Inltlel  hyper- 
ptene  containing  the  four  intermediate  base  points. 

Thus,  with  the  first  estimate  from  (3.4f)  and  the  double 
iteration  process  described  by  equations  (3.4c)  and  (3.4d),  the 
coordinates  of  the  new  field  point  can  be  obtained. 

Step  3 

The  next  step  is  to  determine  the  flow  properties  at  the 
new  field  point.  The  compatibility  equations  (2.15),  (2.16), 
and  (2.18c),  in  finite  difference  form,  are  used  to  determine 
the  flow  properties.  In  order  to  utilize  these  equations,  the 
0-coordinates  and  ^  -coordinate  must  be  determined.  A  rep¬ 
resentation  of  a  portion  of  the  field  point  network  is  depicted 
in  Figure  13.  The  four  intermediate  base  points  are  now  re¬ 
designated  P|^,  P2,  P3>  and  P^  for  simplicity,  while  P^  is  the 
new  field  point.  P^  is  the  point  on  the  initial  hyperplane  from 
which  the  particle  line  through  P^  passes. 

Briefly,  four  flow  properties  at  P^  are  determined  by 
using  equation  (2.15)  written  in  finite  difference  form  along 
the  four  lines  from  the  base  points  to  the  new  field  point. 
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Th«s«  four  •quationt  can  ba  solved  for  v^,  and  p^.  The 
velocity  conponents  give  the  direction  cosines  of  the  general¬ 
ised  velocity  ^  which  is  tangent  to  the  particle  line.  A  line 
with  direction  cosines  equal  to  the  direction  cosines  of  9  is 
then  projected  back  until  it  intersects  the  initial  hyper¬ 
surface.  This  intersection  point  is  P^.  The  compatibility 
equations  (2.16)  and  (2.18c)  along  the  line  in  finite 
difference  form  determine  the  remaining  flow  properties  at  P^. 

The  transformation  needed  to  explicitly  write  the  com¬ 
patibility  equation  (2.15)  has  only  one  restriction  put  upon  it 
by  the  partial  differential  equations  and  it  is  that  the 
characteristic  equation  (2.6)  must  be  satisfied.  It  is  probably 
best  to  take  the  transformation  to  be  linear  and  orthogonal. 

This  does  not  mean  that  the  simplest  form  of  the  required  trans¬ 
formation  is  given  by  the  orghogonal  transformation.  Quite  to 
the  contrary,  nonorthogonal  transfdrmations  could  be  simpler, 
but  they  are  not  easily  determined  analytically  or  numerically. 
For  example,  Thomhill^^  presents  a  nonorthogonal  transformation 
for  two-dimensional  unsteady  flow  which  is  quite  simply  expressed. 
IKis  transformation  was  obtained  by  considering  the  geometry  of 
the  characteristics  surfaces  in  the  three-space.  Unfortunately, 
a  direct  application  of  the  same  procedure  to  the  four- space  is 
not  possible  because  of  the  added  geoeMtrical  complications  of 
the  additional  coordinate. 

Ve  proceed  to  develop  the  linear  orthogonal  transformation 
which  can  be  written  in  matrix  notation  as 
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The  orthogonellty  property  gives 


The  ^^-coordinete  hss  been  taken  normal  to  the  Kach  hyperconoid 

Mfid,  hence,  is  normal  to  the  line  At  the  point  (  i  «  1, 

2,  3,  4).  The  aligned  in  the  direction  along 

16  27 

P^P^.  Fcom  the  work  of  Thornhill  and  Bruhn  and  Haack  ,  it  is 
determined  that  7  0^^  and  7  ^2  should  have  the  following 
forms 

V/fi  =  Css  ^  r  s/n  Pa>%  sift  ^  %in  ^  k 

"  ^ (3.7a) 

VA  =  (|+c»5«Jjr^  ('!' 

^(^+sin<)5i'nt)ik  -f-  ^  T 

where  and  {A  are  parameters,  and  it  is  easily  verified 
by  substitution  that  7  0^  satisfies  the  characteristic  equa¬ 
tion  (2.6).  Also,  it  is  easily  shown  that 


Utilizing  (3.7a-b),  the  complete  transformation  is  determined 
in  the  following  manner.  Equation  (3.5)  is  rewritten  in  the 
following  form 
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l/X*-  =  +  f  D*- 

1/Kz  =  £‘  +  P‘  fS’-+  H‘ 
1/K^  =  L*  +  M‘  +  M‘  +  P‘ 
1/K^‘=  +  S‘  +  T* 


(3.8«) 


(3.8b) 

(3.8e) 

(3.8d) 

(3.8.) 


Thar*  ara  aclll  thraa  dagraaa  of  fraadoa  laft  In  orlantlni 
tba  8*eeordlnataa  I  ao  thraa  tana  can  ba  chosan  arbitrarily. 
Iha  eholeaa  aada  hara  ara 

PxSal  T»0  (3.9) 


Th«  ^  Aligned  with  lo  equmtlng  the 

^2  tens  to  th«  dlroctioa  cosines  of  P^P^  gives 

K.e  s  (Xs-Xi)yz 
KF  ~  (yj-aO/Z 

Ka6  » 

K»H  *■  (Xt~t\)/Z 

vhsrs 

E*=(Xr*)fi)  +(yr*yi) +(*f  (3.10a) 


(3.10a) 

(3.10.b) 

(3.10c) 

(3.10d) 


Noir  utlng  •quatlons  (3.7),  (3*9),  (3*10)  and  the  orthogonality 


conditions ,  the  tarns  in  Q  are  datarmined  with 

aquations 

tha  folloiring 

H  -  V/ai 

(3.1U) 

E-  H 

(3.11b) 

r  -  H  (y»  -  aO/^t'r  -  t'J 

(3. lie) 

(3. lid) 

A  -  E  - 

(3.11.) 

B  -  F  -  ( Vt  /Si) 

(3.11f) 

c  -  s  -  (w./ai) 

(3.11g) 

D  •  -  (Am;  +  BV;  t  Cu<  +  a;  )/v 

(3.11h) 

Q  •  (BS*CF)/(AF-6E; 

(3.111) 

»  -  (CE-A6)/(AF-BEj 

(3.11J) 

L  -  [CF-B6  +(0S-CH)Rj/j5 

(3.11k) 

m-CAS-CE  +(CH-DS)QJ/J5 

(3.111) 

"-[BE-AF  +(0F-BH)<3  +  (AH-PE}R]/J5 

(3.11m) 

where 

ti  *  (BS  'CF)(5  CCE  -/^6)R 


(3.11n) 


These  equations  uniquely  detemine  the  transformation  at  a 
particular  intermediate  base  point.  Using  (3.6)  and  (3.8a), 
it  can  be  seen  that 


^  dij 


§o  ch«t  all  tha  partial  darlvatlvat  froa  tha  tranafoTBatlon  in 
(2*15)  ara  knovn. 


Tha  eoapatlblllty  aquation  (2.15)  la  put  In  flnlta  dlf- 
faranca  fora  by  aubatltutlng  tha  folloirlng  approxlaata  axpras- 
alon  for  tha  partial  darlvatlToa  along  tha  llna  froa  to 


i/S, 


(«r  -  Mi)/K 


(3.12) 


Slallar  aspraaalona  ara  utlllsad  for  tha  darlvatlvat  of  v, 

V,  and  ,;>• 

A  aehaaa  auat  alto  ba  davlaad  to  nuaarlcally  avaluata  tha 
darlvatlvat  In  tha  and  9^-dlractlont.  Thlt  It  ona  of  tha 
fundaaantal  dlffarancat  bataaan  problaat  vlth  tvo  Indapandant 
varlablat  ahara  thara  ara  no  partial  darlvatlvat  laft  to  ba 
avaluatad  at  thlt  point,  and  problaat  with  aora  than  tvo  Indapand¬ 
ant  varlablat*  Thata  partial  darlvatlvat  can  ba  avaluatad  ualng 
tha  following  aquatlona* 


^ 


(3.13«) 


vhara  froa  tha  orthogonality  of  tha  tranaforaatlon 
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(3.13b) 


JX.  X.  ^  ^  ^  ^ 

^  6$m  ^  bA^"^  bh  i An "  6t' 

and  tha  other  partial  darlvatlvat  can  ba  obtained  fron  aquatione 
of  tha  fon 


Aa  a  fir  at  eatiiaata  in  the  iteration  procedure  deacribad  in 
Step  4  balov,  can  he  taken  aa 


l/g'  -  (Ui  +  Mt+  U3  +^4)/^ 


(3.13d) 


In  a  alaiilar  manner ,  the  derivativea  of  v,  w,  and  p»can  be 
dateminad. 

It  ahould  be  pointed  out  that  in  writing  equation  (3.13c) 
it  ia  aaamied  that  the  partial  derivativea  of  the  fora  (^14  (3<  t 

b{A/6\^  t  .  «nd  bu/bt  can  be  approxiaatad  aa 

conatanta  over  the  network.  The  derivativea  bu/ are  not 
conaidared  conatant  over  the  network  in  (3.13a),  however,  bacauae 
tha  tranaforaation  and,  hence,  the  derivativea  given  by  (3.13b) 
are  calculated  for  each  intaraediate  baae  point. 

With  the  equationa  given  above,  all  the  taraa  in  (2.15)  can 
ba  evaluated.  The  last  ten  which  includea  the  chemical  aource 
function  ia  evaluated  at  leaat  initially  at  each  of  the  inter- 
aadiata  baae  pointa  where  all  of  the  flow  propartiea  are  known. 
Thua,  with  i  *  1,  2,  3  and  4,  equation  (2*15)  gives  four  linear 
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algebraic  •quatloaa  which  can  ba  tolvad  for  tha  four  unknowna, 


U5,  Vj,  tad  fy 

Tha  coordinataa  of  tha  baaa  point,  ara  dataminad  by  tha 
iatarsaetion  of  a  lina  through  with  diraction  cosinaa  aqual 
to  thoaa  of  ^  and  tha  hyparplana  through  tha  four  intaraadiata 
baaa  pointa,  P^,  P2>  P3  and  P^.  Tha  aquationa  which  giva  thaaa 
coordinataa,  ara  tha  aquation  of  tha  hyparplana 


+  Pt'  =  1 


(3.14«) 


and  tha  aquationa  of  tha  lina 

— — =.  (3.14b-d) 

U(  Vi  ^ 

Thia  aat  of  linaar  aquationa  ia  aaaily  aolvad  for  tha  coordinataa 
of  fy  Tha  flow  fiald  propartiaa  at  P^  can  ba  dataminad  by 
linaar  interpolation  batwaan  tha  four  intaraadiata  baaa  pointa 

%-A 

in  tha  initial  hyparplana.  Firat,  tranafon  coordinataa  to  tha 

^  coordinataa  whara  tha  t*  coordinataa  ara  all  aqual 
to  tha  aaaa  conatant.  Than,  a  ganaral  flow  proparty  can  ba 
vrittao  in  tha  fom 

^  ~  y*  ^*)  ~  y*  +  Pi  <3.15) 

and  tha  coafficianta  can  ba  dataminad  by  aubatituting  tha 
coordinataa  and  flow  propartiaa  of  tha  four  intaraadiata  baaa 
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points.  Then  the  flow  properties  et  ere  eesily  detenined 
by  substituting  its  coordinates. 

When  the  flow  properties  at  are  known,  the  coapatibility 
equations  along  the  particle  line,  (2.16)  and  (2.18e),  are  used 
to  detemine  the  resuining  flow  properties  at  P^.  For  this  pur¬ 
pose,  these  equations  night  be  put  in  finite  difference  fom 
by  substituting 

Jr"  ~  (3,16.) 

(s,  -  s^)/x 

where 

X  *  -  (»ff  -  ‘‘+  (y  r  (2c  -  Kp  (i'r  -  tl/  (  3*  We) 

(2.16)  and  (2.18c)  can  then  be  solved  for  c^^  and  s^. 

It  should  be  noted  that  in  introducing  the  approxiaations 
of  equations  (3.16),  we  assuae  that  c^  and  s  do  not  vary  nore 
rapidly  iron  the  initial  hyperplanc  to  the  new  point  than  u,  v, 
w,  and  p.  In  many  nonequi librium  flows,  c^  and  s  do  vary  rapidly 
along  the  particle  line  so  that  in  these  cases,  it  nay  be  ad¬ 
visable  to  integrate  equations  (2* 16)  and  (2.18c)  in  several 
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steps  as  proposed  by  Ferri  et.al  and  utilised  in  the  calcula- 

36 

tions  of  Wood  et.al.  The  pressure  can  be  assusMd  to  vary 

linearly  along  the  particle  line  and  the  one -dimensional  or 
streamtube  approximations  allow  a  Runge-Kutta  technique,  for 
example,  to  be  used  to  integrate  the  equations  from  P^  to  P^. 
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Sedney  and  Oerber^^  lave  also  noted  that  utilizing  entropy  in 
nonequl librium  calculations  can  lead  to  large  confutation  errors 
so  that  it  might  be  advisable  in  some  cases  to  use  another 
thermodynamic  variable,  such  as  tenferature,  in  place  of  the 
more  cxistomary  quantity,  entropy.  A  note  of  warning  might  be 
inserted  here  to  indicate  that  the  numerical  Intesratlon  of  the 
reaction  rate  equations,  can  be  difficult  especially  near 
equilibrium  conditions.  Hie  numerical  integration  of  such 
''stiff”  equations  is  Just  beginning  to  bo  handled  adequately 
at  the  present  time. 


Step  4 


The  final  step  is  to  set  \q>  an  iteration  process  to  increase 
the  accuracy  of  the  integration  procedure.  This  is  done  by 
repeating  steps  2  and  3$  but  using  average  values  of  the  flow 
field  properties.  In  equation  (2.13),  for  example,  the  flow 
pTOperties  u^,  v^,  w^,  p^  and  a^  are  replaced  by  average 
values  which  have  the  form 


(ic  J,  (3.17) 


In  equations  (2.l6)  and  (2.l6c),  and  hence,  also  in  (3.14) 
the  following  form  of  average  values  is  used. 


(3.18) 


In  the  iteration  process,  average  values  of  the  flow  properties 
are  used  everywhere  except  in  the  interpolation  procedure  on  the 
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initial  hyptrsurfAc*  and  for  the  initial  conditions  in  the 
integration  of  the  compatibility  equations.  The  iteration 
process  is  continued  until 


(h) 

u,  = 


=  U 


(h-i) 


# 

e 


to  within  prescribed  limits  or  until  they  are  equal  to  the 
ability  of  the  computer  to  determine  them  utilizing  a  finite 
ntaber  of  digits. 

The  iteration  process  described  above,  is  similar  to 

the  modified  Euler's  method  which  is  sosMtimes  referred  to  as 

42 

Heun's  first  method  (see  Lepidus  ,  p.  88)  for  ordinary 
differential  equations.  When  applied  to  ordinary  differen¬ 
tial  equations,  the  truncation  error  is  third-order  in  the 
step  size,  so  that  the  process  includes  terms  of  second  order 
in  the  step  size.  The  process  as  applied  here  to  the  partial 
differential  equations,  does  omit  some  terms  which  are  second 
order  in  the  step  size  through  the  assuoption  that  the  partial 
derivatives  in  (3.13c),  arc  constant  over  the  network.  It 
can  only  be  hoped  that  these  neglected  terms  do  not  appre¬ 
ciably  increase  the  order  of  the  truncation  error  because  it 
appears  that  there  is  no  way  to  include  the  neglected  second- 

order  terms  without  changing  the  network  by  adding  additional 

28 

base  points,  as  Butler  has  done. 

This  completes  the  description  of  the  fourth  step  and 
also  the  entire  field  point  procedure.  Next,  the  procedure 
to  calculate  a  new  point  on  a  body  hyper  surface  is  considered. 
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3.3  The  Body  Point  Procedure 

The  lAM  four  steps  required  for  the  field  point  procedure 
ere  elso  required  for  the  body  point  procedure.  The  body  is 
essuBed  to  be  specified  by  en  equation  of  the  form 

o  (3.19) 

With  this  formulation,  rather  general  body  shapes  and  motions 
can  be  considered.  Ho%#ever,  it  is  assumed  here  that  the  body 
hypersurface  has  no  corners  (i.e.,  no  discontinuities  in  first 
order  derivatives).  Special  techniques  are  required  for 
handling  the  expansion  fans  and  secondary  shock  waves  which 
can  arise  because  of  body  comers. 

Three  base  points  in  the  initial  hypersurface  adjacent 
to  the  body  hypersurface,  and  three  points  which  lie  both  in 
the  body  and  the  initial  hypersurfaces  are  used  to  determine 
a  new  point  on  the  body  hypersurface  in  the  following  steps. 

Step  1 

As  in  the  first  step  in  the  field  point  procedure,  inter¬ 
mediate  base  points  must  be  properly  located,  and  the  flow 
properties  at  them  determined  by  interpolation.  Three  inter¬ 
mediate  base  points  must  be  located  such  that  the  domain  of 
dependence  of  the  difference  scheme  contains  the  domain  of 
dependence  of  the  partial  differential  equations  to  ensure 
ntaerical  stability.  This  was  done  by  inscribing  a  sphere 
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within  the  tetrahedron  In  the  field  point  procedure.  A 
similar  procedure  could  be  proposed  here,  but  it  was  found 
in  doing  two 'dimensional  unsteady  flows,  that  some  experimen* 
tation  was  necessary  in  order  to  properly  locate  the  inter* 
mediate  base  points,  so  that  a  stable  solution  was  obtained. 
Thus,  it  is  probably  better  to  refer  to  the  details  of  the 
t«ro* dimensional  unsteady  procedure  given  in  Chapter  4 ,  Section 
4.3  in  the  description  of  the  PATCH  subroutine,  and  to  state 
that  some  experimentation  might  be  necessary  in  locating  the 
intermediate  base  points. 

Step  2 

The  coordinates  of  the  new  body  point  are  determined  by 
simultaneously  solving  the  equation  of  the  body  hypersurface 
(3*19)  together  with  the  equations  for  three  Mach  h3rpercones 
through  the  three  intermediate  base  points  determined  in  Step 
1.  The  equations  can  be  solved  by  using  the  third*order  Rich¬ 
mond  iteration  procedure  presented  in  Step  2  of  the  field  point 
procedure . 

A  representation  of  a  portion  of  the  body  point  network 
is  shown  in  Figure  14.  ,  P^  and  P^  are  the  intermediate  base 

points,  and  P^,  P^  and  Pg  are  the  base  points  in  the  body  hyper 
surface.  The  new  body  point  is  located  at  ,  while  P^  is  the 
point  of  intersection  of  the  particle  line  from  P|^ ,  projected 
back  to  the  plane  determined  by  P^,  P^  and  Pg. 
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Sfp  3 


Th«  flow  proporciot  at  arc  determined  in  a  manner 
aimilar  to  Step  3  in  the  field  point  procedure.  Three  com¬ 
patibility  equations  can  be  obtained  in  the  form  of  (2.13) » 
eorreaponding  to  the  three  intermediate  base  points.  The 
same  coordinate  transformation  as  used  for  the  field  point 
can  be  used  here.  The  equation  for  the  boundary  condition 


4  11^7  4  =  o 

6^  dl  6t' 


(3.20) 


replaces  the  fourth  compatibility  equation  used  in  the  field 
point  procedure.  One  slight  change  is  necessary  to  evaluate 
the  derivatives  in  the  0^-  and  0j^-dircctions  in  (2.13) 
because  only  three  and  not  four  intermediate  base  points  are 
available.  One  of  the  base  points  on  the  body  h3rpersurface 
must  be  used  as  the  fourth  point.  Prom  the  three  coeq>at- 
ability  equations  and  the  boundary  condition,  the  flow 
i'prpperties  u^,  V|^,  Wj^  and  p^^  arc  determined. 

Meat ,  the  particle  line  on  the  body  hypersurfacc  must  be 
projected  back  to  obtain  its  intersection  with  the  initial 
hypersurface.  The  three  body  base  points,  P^,  P^  and  Pg  in 
Figure  14,  determine  a  plane  (not  a  hypcrplanc)  which  is  ex¬ 
pressed  in  the  following  form. 


Ax  Cy  +  Ds  -  1 
Ey  +  Fx  +  Gt'  -  1 


(3.21a) 

(3.21b) 
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The  fix  coefflcicnCe  cen  be  determined  by  eubstitutlng  the 
coordinetes  of  the  three  points  into  the  equations,  end 
solving  the  resulting  set  of  six  equations.  This  plane  is 
an  approximation  to  the  surface  which  is  the  intersection  of 
the  initial  hypersurface  and  the  body  hypersurfacc.  The  point 
P-,  which  is  the  intersection  of  the  particle  line  with  the 
plane  (3.21),  is  determined  by  considering  a  second  plane  con¬ 
taining  P^,  the  generalized  velocity  ^  and  the  vector  normal 
to  the  body  hypersurface  at  P^^ ,  %#hich  is  given  by 


N 


t 


6i  dt 


(3.22) 


This  plane  can  also  be  expressed  in  the  form 

Hx  +  ly  -»■  Jz  •  1  (3.23s) 

Ky  +  Lz  :  Nt‘  -  1  (3.23b) 

The  coordinates  of  P^  are  obtained  by  solving  equations  (3.21) 
and  (3.23).  The  intersection  of  the  two  planes  results  in  a 
unique  point,  because  the  plane  (3.23)  contains  N|^,  which  is 
normal  to  the  body  hypersurface.  Hence,  the  intersection  of 
(3.23)  with  the  body  hypersurface  is  the  particle  line,  and 
because  (3.21)  also  lies  within  the  body  hyper sur face ,  the 
intersection  is  a  single  point. 

After  the  particle  line  base  point  P^  is  determined,  the 
flow  properties  at  it  are  determined  by  interpolation  in  the 
initial  hypersurface.  The  compatibility  equations  (2.16)  and 
(2.18c)  along  the  particle  line  are  used  to  determine  the 
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rMuinlng  flow  proportiot  at  tha  now  body  point  in  tht  taaa 
■ancar  aa  in  tha  fiald  point  procadura. 

Stap.  » 

Stapa  2  and  3  abova  ara  itaratad  to  convarganca  in 
axaetly  tha  aaaa  nannar  aa  in  tha  fiald  point  procadura. 

3.4  Tha  Shock  Point  Procadura 

Tha  ahock  point  procadura  ia  tha  noat  coaplicatad  of  tha 
thraa  procaduraa  conaidarad  hara.  Briafly,  tha  ahock  point 
procadura  can  ba  auHurizad  aa  followa.  Givan  a  point  which 
liaa  both  on  tha  ahock  and  initial  hyparaorfacaa ,  togathar  with 
tha  eoaplata  flow  fiald  incidant  on  tha  ahock,  tha  location  of 
a  new  point  on  tha  ahock  hyparaurfaca  wuat  ba  found.  Tha  flow 
proper tiaa  behind  tha  ahock  at  tha  new  point  and  tha  orien¬ 
tation  of  tha  ahock  hyparaurfaca  through  tha  point  uuat  alao 
ba  datarsinad.  Aa  in  tha  body  point  procadura,  thraa  baaa 
polata  in  tha  initial  hyparaurfaca  adjacent  to  tha  ahock 
hyparaurfaca  ara  aaauaad  givan.  Actually  only  one  ahock  baaa 
point  which  liaa  both  on  tha  ahock  and  initial  hyparaurfacaa  ia 
required  to  datamina  tha  new  ahock  point,  but  an  additional  two 
ahock  baaa  pointa  will  ba  required  in  Stap  1  to  datamina  inter - 
■adlata  baaa  pointa  aa  waa  nacaaaary  for  tha  body  point. 

Tha  ahock  point  procadura  haa  tha  aaaa  four  atapa  aa  tha 
fiald  and  body  procaduraa:  (1)  locate  and  datamina  tha  flow 
propartiaa  at  tha  intamadiata  baaa  pointa,  (2)  locate  tha 
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position  of  the  new  point,  (3)  celculetc  the  flow  properties 
St  the  point,  end  (4)  iterate  the  solution  to  second  order 
in  the  step  size. 

Step  I 

This  step  is  almost  exactly  the  sasie  as  Step  1  of  the 
body  point  procedure.  Also,  as  in  the  body  point  solution,  a 
detailed  description  of  a  process  to  locate  the  intermediate 
base  points  could  be  given  here,  but  because  some  experimen¬ 
tation  will  probably  be  necessary  to  assure  numerical  stability, 
the  reader  is  referred  to  procedures  used  for  two-dimensional 
unsteady  flow  in  Chapter  4,  which  probably  can  be  sinf^iy 
generalized  to  the  four  variable  case. 

Stgp  2 

A  representation  of  a  portion  of  the  shock  point  network 
is  shown  in  Figure  I5.  ?2  and  are  the  three  inter¬ 
mediate  base  points  %rhich  were  determined  in  Step  1.  P^^  is 

the  base  point  which  lies  in  both  the  shock  and  initial  hyper- 
surfaces.  The  new  shock  point  that  is  to  be  determined  is 
labeled  P  . 

The  first  estimate  for  the  location  of  P^  is  obtained  by 
determining  the  intersection  of  the  shock  hypersurface  and  the 
three  Mach  hypercones  from  the  intermediate  base  points.  The 
shock  hypersurface  is  initially  approximated  by  the  hyperplane 
normal  to  and  passing  through  P^.  is  the  unit  vector 
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Th«  shock  hypsrpisns  Chen  has  the  form 


+  +  0  (3.25) 

The  coordinates  of  can  be  obtained  by  using  the  third  order 
Richmond  iteration  scheme  as  in  the  field  and  body  point 
procedures . 

St«p  3 

The  procedure  for  determining  the  flow  properties  behind 
the  shock  hypersurface  at  the  new  shock  point  is  more  involved 
then  the  corresponding  procedures  for  the  field  and  body  points 
First,  notice  that  the  shock  hypersurface  at  has  been  approx 

imeted  by  the  hypcrplane,  (3.25) >  >o  that  in  the  first  step  of 
the  interation 

The  flow  incident  on  the  shock  wave  at  P^  is  given,  so 

V54..  ’'5+.  P5+  Cn5+  •re  known  or  can  be  deter¬ 

mined  by  interpolation.  The  plus-sign  subscript  indicates  a 
flow  property  in  the  incident  stream  while  a  minus-sign  sub¬ 
script  iadieates  flow  properties  just  behind  the  shock  hyper- 
smrface.  The  shock  wave  equations  (Rankine-Hugoniot  equations) 
are  need  to  determine  the  properties  behind  the  shock. 


Introducing  the  following  notation 


“  (Ut  ^  ^  ^  ^  ^v) 

where 

)/(^V  ■*■  ^s*) 

and 

"  H J ^ ^  \ 

the  aquations  for  a  noving  shock  wave  can  ba  written  in  the  fora 


(3.26d) 

f=v  +  ^ 

(3. 26.) 

(3.26f) 

'  ^T. 

(3.26g) 

Ch+  =  (^n  =  1 ,  .  .  .  , 

(3.26h) 

State  equation  is  also  needed. 

h  =  H  (f>,  /),  Cl, - ,  Cm) 

(3.261) 

These  equations  cannot  be  solved  explicitly  for  the  properties 
behind  the  shock  in  terms  of  the  properties  in  front  of  the 
shock,  because  the  state  equation  cannot  always  be  written 
explicitly  in  the  form  of  an  equation.  The  shock  wave  equations 


(3.26a) 

(3.26b) 

(3.26c) 
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■tttt  b«  solved  nvMrically  for  a  general  problsa,  such  as 

44 

tha  flow  of  a  real  gas  (a.g.,  for  air,  saa  Moackal  ). 

After  q||,  has  bean  datarwlnad  using  (3*26),  tha  velo¬ 
city  casq>onants  are  given  by  tha  following  aquations. 

M.  -  W*  +  (3.27«) 

v:  -  v;+  N.JJ  $  (3.27b) 

KT.  -  \^4  ^  (3.27c) 

Thua,  all  the  flow  properties  at  behind  tha  shock  hyper- 
■orfacas  can  be  dataminad. 

Step  4 

Equations  (3.26)  and  (3*27)  datanina  an  astinata  of  tha 
flow  properties  at  P.  behind  tha  shock  wave,  biit  there  are 
three  coeipatibility  equations  which  hold  along  the  lines 
Pj^P^,  P2P^  and  In  general,  tha  properties  dataminad 

above  will  not  satisfy  these  cowpabtlity  aquations.  However, 
a  new  astinata  for  tha  shock  hyparsurfaca  nomal  vector,  N 
can  be  nada  such  that  tha  conditions  behind  the  shock  tend 
toward  tha  satisfaction  of  tha  conpatibility  aquations. 

Tha  new  astiwata  of  11.^  I*  obtained  as  follows.  Tha  flow 
properties  which  appear  in  tha  cowpatibility  aquation  (2.13) 
are  expanded  in  Taylor  series, which  is  truncated  after  tha 
first  ordsr  derivatives 
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(3.28) 


where  ■inue-elgn  eubecripte  heve  been  onitted.  Siniler  eque- 


Clone  can  be  written  for  v. 


(n+l)  (n+l)  (n+1) 

»  5  5 


eetlwete  of  le  given  by 

*•  •  J  ^  \  J 


(3.29) 


Tta*  thra*  corraction  Cana  ,  /Olitl 

■*!>  •Jj 


end  AM: 


ere 


determined  by  eubetituting  (3.28)  end  ite  connterperte  for 
V,  w,  end  p  into  the  three  competlbillcy  cquetione,  (2.15)» 
correeponding  to  •  ?2  end  P^.  The  pertiel  dcrivetivee  of 
u,  V,  w,  end  p  with  reepect  to  N  M.^  end  N..  auet  be 
determined  before  the  correction  terme  cen  be  obteined. 

Theee  derivetivee  aey  be  obteined  from  equetione  (3.26)  end 
(3*27),  but  they  cennot  be  determined  enelyticelly  unleee 
the  etete  equetion  cen  be  written  explicitly.  They  cen  be 
determined  nijmericelly ,  end  perhepe ,  this  if  the  simpleet  wey 
even  when  they  cen  be  determined  enelyticelly.  For  exemple, 
e  derivetive  could  be  celculeted  using  the  following  epproxime* 


tion. 


where  ^  *  emell  ohenge  in 


After  the  three 


correctionu  ere  determined,  the  fourth  correction  is  obteined 
from 


65 


(3.30) 


(Ns»j.+ +  (Ns^  +  ^W»3s-) 


This  proctdurt  is  rcpastsd  in  sn  Itsrstion  process  ss 
la  ths  field  end  body  point  procedures.  In  Step  I,  is 
replaced  by  sn  sversge  value  given  by 


C=iK*0"^vOj 
+( V  ^  O'] 


(3.31) 


uhere  is  given  by  (3.29).  In  Step  2,  note  that 

is  not  equal  to  S  .  given  by  (3.29).  Average 

*0 

values  are  used  in  the  hypercone  and  coupatibility  equations 
Just  as  in  the  field  and  body  point  procedures. 


3.*^  The  Initial  Value  Hypersurface 

The  Cauchy  problem  being  considered  here  requires  thet  the 
flow  properties  be  completely  specified  on  an  initial  hyper¬ 
surface.  This  means  that  a  three-dimensional  flow  field  must 
be  specified  before  the  complete  four -dimensional  flow  field 
can  be  solved.  By  no  means,  is  this  a  trivial  problem.  Of 
^course,  the  type  of  initial  value  hypersurface  required  for 
a  given  problem  will  depend  on  the  specific  problem.  A  few 
of  the  properties  which  seem  unique  to  the  general  method  of 
characteristics  and  which  might  be  used  to  determine  the 
initial  hypersurface  are  mentioned  in  the  following  discussion. 
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In  most  practical  problams  of  unsteady  three "dlBcns tonal 
flow,  the  initial  hypersurface  is  the  initial  steady  three- 
dieensional  flow.  Thus,  the  steady  flow  is  required  before 
the  unsteady  flow  can  be  determined.  This  brings  to  mind, 
specializing  the  general  method  of  characteristics  to  solve 
the  three-dimensional  steady  flow.  This  is  discussed  in  the 
next  section.  Only  steady  supersonic  flow  can  be  solved  by 
the  method  of  characteristics,  because  the  equations  of  steady 
subsonic  flow  are  not  hyperbolic,  but  this  is  one  way  in  which 
the  initial  hypersurface  might  be  determined  for  supersonic 
flow. 

Certain  steady  flows  might  also  be  solved  in  the  follow¬ 
ing  manner.  Consider  the  hypersonic  flow  about  a  blunt  nosed 
body  as  depicted  in  Figure  16.  The  flow  in  the  initial  hyper¬ 
surface,  t'  «  0,  Buiy  be  only  an  estimate  of  the  exact  flow.  The 
flow  can  then  be  calculated  in  the  t* -direction  with  the  body 
surface  unchanged  with  time.  The  initial  estimate  oan  be  con¬ 
sidered  physically  as  the  result  of  soaM  body  motion  for  t'  <  0. 
After  a  certain  amount  of  time  has  passed,  the  transient  decays 
and  steady  flow  is  attained  at  t'  *  t*  .  For  t'  >  t'  ,  the 
flow  field  does  not  change  with  time.  In  this  way,  the  steady 
flow  can  be  calculated  starting  from  only  an  estimate  of  the 
flow  field.  Note  that  the  initial  estimate  must  not  violate 
the  equations  of  change  of  the  gas.  Hence,  the  initial  esti¬ 
mate  is  not  arbitrary,  but  must  be  physically  reslizable. 
Obtaining  such  an  estimate  for  nonequilibriixD  flow  may  not  be 
simple,  but  perhaps  rather  crude  Initial  estimates  would  still 
yield  the  physically  correct  results.  This  could  be 
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The  difficulty  in  obtaining  an  initial  catiaatc  which  ia 
coaipatibla  with  tha  eqiiationi  of  change  could  be  avoided  by 
uaing  the  following  technique.  Start  with  an  initial  steady 
flow  which  ia  kaown  exactly.  For  exaaple,  the  hypersonic  flow 
about  a  spherically  capped  cylinder  wight  be  taken  as  the  known 
steady  flow  as  shown  in  Figure  17.  The  spherical  cap  can  then 
be  distorted  into  som  other  shape,  perhaps  ellipsoidal,  and 
after  the  transient  decays,  the  steady  flow  about  the  dis¬ 
torted  body  is  obtained.  In  the  sawe  way,  the  known  initial 
steady  flow  slight  be  taken  as  the  supersonic  flow  about  a  cone 
at  sero  angle  of  attack.  The  cone  can  then  be  rotated  to  a 
finite  angle  of  attack,  and  after  the  transient  decays,  the 
steady  flow  about  a  cone  at  an  arbitrary  angle  of  attack  is 
obtained.  In  this  way,  the  steady  supersonic  flow  about  sowe 
rather  arbitrary  bodies  can  be  calculated. 

It  should  be  sMntioned  that  the  general  wethod  of 
characteristics  for  unsteady  flow  applies  to  subsonic  as  well 
as  supersonic  flows.  Given  the  initial  steady  subsonic  flow, 
the  flisw-field  history  can  be  cowputed  within  the  doeMin  of 
deteminacy  of  the  initial  data.  The  dowain  of  deterwinacy  is 
lisiited  in  the  t' -direction  by  the  awoxmt  of  initial  data  given 
in  the  initial  hypersurface,  but  the  calculation  can  be  carried 
as  far  in  the  t* -direction  as  desired,  nerely  by  enlarging  the 
region  of  the  initial  data.  For  practical  reasons  discussed 
in  Chapter  3,  the  wethod  of  characteristics  for  unsteady  flows 
will  probably  never  be  used  for  subsonic  flows  unless  the  flow 
is  desired  for  only  a  very  short  period  of  tiae. 
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3.6  Siapllfylng  th«  CenTtl  Proc*dttre 


It  will  not  Always  b«  nacassary  to  consldar  tha  thraa- 
dimensional  unstaady  flow  of  a  nonaquilibriun  gas.  Tha  wethod 
of  charactaristics  can  also  ba  appliad  to  siaplar  problaaa. 
Heans  for  spacia Using  tha  gadaral  procadure  can  ba  vary  use¬ 
ful  and  should  ba  kept  in  Bind  as  tha  general  procedure  is 
developed.  Specifically,  three  simplifications  are  considered 
hare.  They  are  listed  below,  together  with  a  brief  discussion 
of  how  tha  siaplification  is  carried  out. 

Nonequilibrium  flow  simplified  to  equilibrium  flow. 

In  this  case,  the  simplification  is  accomplished  by  eliminating 
the  following  variables  and  operators  :  DCj^/Dt' 

and  d  Also,  the  equilibrium  spaed  of  sound,  a^ ,  should 

be  substituted  for  the  frozen  speed  of  sound.  This  should  be  a 
rather  simple  process. 

Three -dimensional  unsteady  flow  simplified  to  two- 
dimensional  unsteady  flow.  Here,  one  of  the  space  coordinates 
is  to  be  dropped  so,  for  instance,  the  following  terms  are  sat 
equal  to  zero:  z,  d  /dz  ,  w,  3^^  and  ^  In  the  co¬ 
ordinate  transformation,  is  also  set  equal  to  zero.  This 

too  seams  to  ba  a  rather  simple  process. 

Three-dimensional  unstaady  flow  simplified  to  three- 
dimensional  steady  flow.  In  this  case,  the  specialization  is 
not  easily  accomplished.  First,  time  must  be  eliminated  from 
the  problem  by  setting  t'  and  d  /dt'  equal  to  zero.  The 
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character Ittic  aquation,  (2.6)  changes  to 


(3.32«) 


•0  that 


M*’  vr*’-f  w- 


(3.32b) 


for  the  problea  to  remain  hyperbolic.  The  equation  of  the 
Mach  cone  changes  to 

+  (m^+ V*--  v-Cif'if.'jCy- aO  (3.32c) 

Vote  that  (3.32c)  is  not  obtained  from  (2.9)  merely  by  setting 
t'  equal  to  aero.  In  order  for  the  Mach  hypercone  to  be  written 
in  a  general  form  so  that  it  can  be  easily  simplified,  the 
characteristic  equation,  (2. 6)^  must  be  written  in  the  form 

i  i  k.  ^  -  o 

n.i  M»1  6i/,„  do<„ 


(3.32d) 


where 


and  then  the  Mach  hypercone  can  be  mritten  in  the  form 

Z  t  “fh;)  =  o 

h»t 


(3.32e) 


where  is  the  matrix  made  up  of  the  cofactors  of  the  deter¬ 
minant  of  b  .  When  t'  is  set  equal  to  sero  in  (3.32d),  equation 
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(3*32«)  will  give  eh*  correct  equation  for  the  Mech  cone. 

Another  coeplicetlon  in  epecielizing  to  eteedy  flow  ie  thet 
the  coordinate  trenefometion ,  (3.8e),  cannot  be  easily  siai- 
plified  to  steady  flow.  Hence »  a  different  traneformation  such 
as  the  one  used  by  Ferri^^  and  Fowell^^  is  required  for  the 
steady  flow.  Because  of  the  above  difficulties,  it  night  be 
easier  to  develop  a  procedure  specifically  for  three -dinensional 
steady  flow,  such  as  the  %rork  reported  by  Fowell,  than  to  try 
to  specialize  three -diMnsional  unsteady  flow.  This  specula¬ 
tion  can  be  verified  when  the  general  procedure  is  prograaned 
for  the  digital  computer. 
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CHAPTER  4 

THEPROGRAMMIHG  OF  THE  PROCEDURE 


In  thit  chapter,  tha  datalla  of  the  progranning  of  the 
Mthod  of  charactarittlca  are  praaantad.  Thla  Ih  the  portion 
of  thla  study  which  pertaina  to  the  second  part  of  the  object 
stated  initially,  that  is,  the  proving  (or  disproving)  of  the 
feasibility  of  the  eethod  by  prograesiing  end  utilizing  the 
Mthod  to  solve  practical  flow  fielda.  This  is  by  no  Mans  a 
■inor  portion  of  the  total  project. 

4,1  Practical  Considerations 

Many  factors  were  considered  in  deciding  exactly  what 
programs  could  be  feasibly  attempted  with  presently  available 
high  speed  electronic  digital  computers.  First,  the  limitations 
of  the  IBM  7094  to  be  utilized  for  the  calculations  had  to  be 
considered.  This  computer  has  approxUiately  32,000  words  of 
Mgnetic  core  storage  available.  About  half  of  this  storage 
should  be  set  aside  for  the  programs,  which  leaves  17,000  words 
for  the  storage  of  data  on  the  initial  hypersurface.  It  has 
been  found  that  the  new  data  can  be  written  over  the  initial 
data  as  the  calculation  progresses  if  a  few  temporary  storage 
locations  are  used.  In  this  wsy,  no  more  than  one  hypersurface 
of  data  needs  to  be  stored  at  one  tiM.  A  few  simple  calcula* 
tions  can  check  the  adequacy  of  this  storage  for  various 
problems.  For  three -diMns ions  1  unsteady  nonequi librium  flow, 

(9  -f  ■)  words  of  data  (the  values  ofx,  y,  z,  t',  u,  v,  w,  p, 
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a,  «nd  c^)  Buat  b«  ator«d  for  each  point.  If  N  la  tight,  for 
tztaiplt,  data  could  ba  atorad  for  1000  pointa.  In  the  initial 
hyparaurfact ,  thia  would  give  10  pointa  in  a  typical  linear 
diaanaion,  which  hardly  aaema  adeqxiate  to  deacriba  a  practical 
flow  field.  If  equilibrium  flow  ia  conaidared,  1900  pointa 
could  be  stored  for  the  initial  hyperaurface ,  but  atill  thia 
would  give  only  12  pointa  in  a  typical  linear  diaanaion.  Finally, 
conaider  a  two*dimenaional  unateady  equilibrim  flow  which  re- 
quiraa  aeven  words  of  data  at  a  point.  This  would  allow  23OO 
points  in  the  initial  surface,  or  5O  points  in  a  typical  linear 
dimension.  This  ia  a  more  reasonable  number  of  points  to  describe 
a  flow  field. 

Another  machine  limitation  ia  the  operating  speed  of  the 
computer.  Even  though  calculations  can  now  ba  carried  out  in  a 
time  on  the  order  of  microseconds,  programaMra  have  been  able  to 
demand  more  and  more  of  computers , as  in  the  problem  considered 
here.  Computation  time  must  be  considered,  because  this  is  what 
determined  the  financial  cost  of  carrying  out  the  calculation. 

A  quantitative  estimate  of  operating  time  is  not  possible  before 
the  programs  are  written,  but  it  is  possible  to  draw  some  quali¬ 
tative  conclusions.  One  result  from  two-dimensional  steady  flow 
calculations,  is  that  nonequilibrium  and  real  gas  thermodynamic 
subprograms  tend  to  require  relatively  larger  amounts  of  computer 
operating  time.  The  complicated  iterations  discussed  in  Chapter 
3  will  require  large  amounts  of  operating  time,  but  these  cannot 
be  eliminated  without  reducing  the  calculations  to  rather  trivial 
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probiMt.  Thus,  it  sppssrt  advisable  Co  consider  a  perfect  gas 
as  a  first  case  Co  keep  operating  time  reasonable. 

The  programs  were  mritten  in  Che  FORTRAN  II  computer 
language  by  Mark  Sustman,  and  the  author.  Neither  of  the  pro- 
grasMera  was  familiar  with  the  FORTRAN  language  when  Che  project 
was  begun.  SoaM,  or  perhaps,  most  of  Che  programs  if  mriccen  in 
a  more  machine  oriented  language ,  such  as  FAP  would  require  less 
storage  and  less  operating  time.  However,  due  to  the  inexperience 
of  the  programsMrs,  it  was  decided  that  the  more  mathematically 
oriented  FORTRAN  language  could  be  utilised  in  a  much  shorter 
time  period. 

Vith  the  above  sMntioned  considerations  in  mind,  it  was 
decided  that  the  general  problem  formulated  in  Chapter  3  must 
be  limited  slightly,  in  order  to  write  sosm  useful  programs  at 
the  prestnt  time.  These  limitations  have  been  incorporated  in 
such  a  way  that  the  difficulty  in  removing  them  should  not  be 
great ,  or  should  not  require  any  fundamental  changes  in  the 
programs.  First,  the  problem  was  limited  to  two-dimensional 
unsteady  flow  as  discussed  in  Section  3*6.  This  reduces  the 
storage  required  to  a  reasonable  quantity,  and  is  rather 
simply  carried  through.  This  restriction  is  also  rather  easily 
rsmoved  and  in  most  equations  (and  therefore,  also  in  program 
statesMBts),  the  neglected  terms  which  would  contain  the  third 
space  coordinate  can  be  re-introduced  by  Inspection. 

To  keep  computing  times  reasonable,  it  was  decided  to 
consider  only  a  perfect  gas  in  this  initial  study.  Equilibrium 
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r«al  gat  CharaodynaBlct  can  ba  Introduced  Mraly  by  replacing 
the  perfect  gat  thexmodynainic  tubroutinet  with  their  real  gat 
counttrpartt .  The  addition  of  nontquilibriue  ef facta  vill  re¬ 
quire  different  cheuical  and  tharaodynaaic  prograat,  and  the 
■odification  of  touM  of  the  prograat  to  include  the  integration 
of  the  tpecidt  continuity  equation  along  the  ttreaelioe. 

There  it  alto  an  argueent  for  tiaplifying  the  probleu, 
other  than  becaute  of  progranming  reatrictiont .  The  problem 
thould  be  aiuplified  at  much  at  possible  while  still  retaining 
its  major  features,  so  that  complexities  from  auxiliary  sources 
such  at  thermodynamic  tubroutinet  do  not  overshadow  the  more 
important  effects  which  are  to  be  investigated.  Thus,  initially 
one  space  dimension  can  be  dropped  and  a  perfect  gas  considered. 
Then,  after  the  programs  for  this  problem  are  working  the  further 
complexities  of  the  third  space  dimension  and  more  complicated 
thermodynamics  can  be  added. 

The  problem  to  be  considered  is  still  very  challenging, 
becaute  three  independent  variables  are  to  be  considered.  Few 
people  as  yet  have  had  success  in  using  the  method  of  character- 

28 

istics  for  problems  in  three  independent  variables.  Only  Butler 
and  Talbot^ utilizing  Butler's  formulation  have  considered  two- 
disMnsional  unsteady  flow.  Butler's  approach  is  not  exactly  the 
method  of  characteristics  as  pointed  out  in  Section  3*1>  and  it 
is  not  obvious  that  it  can  be  easily  extended  to  four  independent 
variables.  The  resulting  programs  will  also  allow  the  checking 
of  some  of  the  proposed  calculation  schemes  discussed  in  Section 

3.5. 
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4.2  Th«  G€n«ral  Approach 

Thr««  tubroutin*!  th«t  datcrmin*  th«  Ioc«tion  and  flow 
proper tl«f  at  the  field,  body, and  ahock  points  are  the  heart 
of  the  prograaned  procedure.  These  are  the  most  important  sub¬ 
routines  and  all  of  the  other  subroutines  and  the  executive  pro¬ 
gram  (or  main  program,  as  it  is  termed  in  programing  manuals) 
merely  support  and  aasiat  these  subroutines  in  calculating  the 
flow  fields. 

The  executive  program  organizes  the  points  in  the  initial 
surface  and  calls  the  three  principle  subroutines  to  calculate 
new  points  on  the  next  surface.  It  then  accepts  the  results 
from  the  subroutines  and  organizes  the  data  in  the  next  surface. 
These  cucfaces  will  be  referred  to  as  time  surfaces,  because 
they  are  usually  close  to  being  t*  •  constant  planes.  The  more 
general  a  type  of  flow  the  executive  program  can  handle;  the 
more  complicated  and  sophisticated  it  must  be.  In  this  Initial 
study,  the  executive  was  kept  as  simple  as  possible  and  because 
of  this,  it  is  very  much  oriented  to  the  particular  type  of  flow 
field  being  solved.  Hence,  if  a  completely  different  type  of 
flow  is  to  be  calculated,  a  new  executive  program  i#ould  be 
required . 

It  should  be  pointed  out  that  the  programs  were  written 
with  the  three  basic  subroutines  for  field,  shock,  and  body 
points  being  supported  by  all  other  programs,  because  these 
are  the  fundamental  elements  In  the  calculation,  and  are 
universal  to  all  Invlscld  flow  problems.  Thus,  they  can  be 
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utlliztd  in  All  other  flow  cAlculatlonf  of  Chit  type,  aerely 
by  eupp lying  the  proper  eupportlng  progremf.  If  the  celculetion 
of  A  new  field,  ehock,  end  body  pointe  hed  been  combined  in 
lerger  progreme  together  with  other  Aupporting  CAlculetione , 
oriented  towerd  specific  flow  fields  end  boundsry  conditions, 
they  could  not  heve  been  utilized  as  reAdily  in  different 
flow  cAlculAtions. 

The  general  progran  organization  is  shown  in  Figure  18, 
in  the  fom  of  a  calling  sequence  diagram.  The  three  sub¬ 
routines  to  calculatie  field,  body,  and  shock  points  hare  been 
given  the  mnemonic  names  FLDPT,  BDYPT,  and  SHKPT.  The  programs 
can  be  divided  into  four  general  sections  as  shown  in  the 
figure.  The  executive  is  the  first  section.  In  the  second 
section,  there  are  four  subroutines  which  control  interpolation, 
the  number  of  points  in  a  surface ,  and  reading  or  writing  data 
on  the  tape.  The  third  section  has  the  subroutines  which  actually 
calculate  new  data  points.  In  the  last  section,  there  is  only 
one  subroutine  which  does  an  auxiliary  interpolation.  Three 
FORTRAII  functions  are  also  utilized.  One  is  named  THETA,  which 
is  called  by  TRUNET  and  PATCH.  The  last  two  functions  carry 
out  the  thenBod3mamic  calculations  and  are  called  ADET  and 
ROEOET.  A  more  detailed  discussion  of  each  of  the  subroutines 
is  given  in  the  next  section.  ^  * 
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Data! led  Dlicufloni  of  the  Progr— ii 

In  the  following  discustiona ,  only  those  techniques  and 
procedures  which  have  not  been  discussed  in  Chapter  3  will  be 
presented.  General  flow  diagrams  and  complete  program  listings 
ere  presented  in  Appendix  C.  These  discussions  are  intended 
to  point  out  the  function  of  the  programs,  tho  checking  pro¬ 
cedures  built  into  them  and  any  unusual  or  special  techniques 
which  were  found  necessary  to  insure  successful  operation.  A 
knowledge  of  the  FORTRAN  II  programning  language  and  the  IBM 
7090/94  Data  Processing  Systems  is  assumed  (see  References  46 
and  47). 

Executive  Program 

As  mentioned  in  Section  4.2,  the  executive  program  is 
oriented  toward  solving  a  particular  flow  field.  The  executive 
considered  here  and  listed  in  Appendix  C  is  for  the  flow  in 
the  subsonic,  transonic  and  supersonic  regions  between  a  de¬ 
tached  shock  wave  and  the  surface  of  a  blunt  body  moving  at 
high  supersonic  speed.  This  program  is  Xn  improved  version  of 
an  earlier  program  which  was  used  to  calculate  the  flow  over 
a  wedge  in  .i  s\ipersonic  or  hypersonic  flow. 

This  executive  reads  the  data  for  the  initial  surface 
from  a  tape  where  it  had  been  previously  written  in  the  proper 
order  and  format.  The  data  is  stored  in  the  common  portion  of 
the  magnetic  core  storage  in  three  large  arrays,  named  FLPTSl, 
BDFTSl  and  SKPTSl  for  field,  body,  and  shock  points,  respectively. 
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The  executive  merely  treneferi  variouf  wor^e  of  data  from  theae 
locatlona  to  other  locations  in  the  common  storage,  where  the 
PATCH,  FLDPT,  BDYPT  and  SKKPT  subroutines  pick  xip  the  input  data 
for  their  calculations. 

Four  data  cards  are  read  by  the  executive  to  specify 
auxiliary  parameters.  The  first  card  contains  niabers  which 
specify  what  the  program  is  to  do  on  a  particular  run,  as  ex¬ 
plained  in  the  listing.  Free  stream  conditions  are  specified 
on  the  second  card  and  the  nondimensionalizing  velocity,  U, 
together  with  the  convergence  test  values  are  specified  on  the 
third  and  fourth  cerds. 

The  organization  in  the  executive  calculates  data  on  an 
odd  nimibered  surface  which  contains  only  field  points.  It  then 
uses  these  points  together  with  the  shock  and  body  points  from 
the  initial  surface  to  calculate  field,  shock,  and  body  points 
on  an  even  nunbered  surface.  The  detailed  organization  of  the 
field  points  is  shown  in  Figure  19.  The  new  shock  and  body 
points  are  calculated  using  two  field  points  from  the  rows 
closest  to  the  shock  and  body  surface  on  the  odd  nissbered  surface. 
The  body  and  shock  points  required  as  base  points  are  taken  from 
the  even  numbered  surface. 

After  the  executive  has  sequencially  calculated  the 
specified  nunber  of  new  time  surfaces,  it  %n:ites  the  final  data 
on  tape.  This  data  is  in  such  a  form,  that  the  executive  can 
read  it  from  the  tape  at  a  later  time  and  continue  the  calculation. 
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SCAN 


Thif  is  a  FAP  coded  subroutine  i#hich  aerely  positions  the 
dete  tape  so  that  the  proper  data  can  be  read  froei  it  or  written 
upon  it.  The  data  is  written  on  the  tape  in  logical  records 
separated  by  end  of  file  marks.  Each  logical  record  is  the 
data  for  one  initial  surface.  SCAN  has  one  integer  argument 
which  specifies  the  logical  record  on  the  tape  which  is  to  be 
read  or  written. 

ADDROW 

As  can  be  seen  from  Figure  19,  one  row  of  points  is  lost 
for  each  step  taken  to  a  new  surface.  It  is  possible  that  in 
certain  flows ,  such  as  the  blunt  body  case  considered  here , 
that  these  points  will  be  spread  over  a  larger  area  at  each 
step.  The  density  of  the  points  will  then  decrease.  ADDROW 
is  used  to  add  more  rows  of  points  between  the  existing  boundary 
rows  by  linear  interpolation.  The  rows  must  have  been  evenly 
spaced  by  TRUNET  before  ADDROW  is  called,  because  ADDROW  assiaes 
even  spacing  in  doing  its  calculation , though  slignt  variation  in 
spacing  can  be  tolerated  by  ADDROW. 

TRPNET 

The  location  of  points  on  a  new  surface  is  a  function  of 
the  spacing  of  the  points  on  the  initial  surface,  the  velocity 
field  and  the  speed  of  sound  throughout  the  field.  Because  the 
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late  two  factor*  cannot  be  controlled,  the  locationa  of  the 
point!  in  the  new  surfacea  tend  to  drift  out  of  their  initial 
configuration  after  a  few  atepa  to  new  tiae  aurfacea  have  been 
made.  TRUNET  uaea  linear  interpolation  to  "true  up"  the  net 
ia  both  tiae  and  apace.  It  can  be  called  periodically  after  a 
few  atepa  to  reatore  the  net  to  ita  original  configuration. 

TRUNET  ia  juat  a  alightly  different  veraion  of  the 
executive  program.  It  uaea  an  organization  aimilar  to  the 
executives  which  ia  ahown  in  Figure  19,  but  it  calculates  new 
points  in  a  slightly  different  sequence,  such  that  it  can  inter 
polate  bet%#een  two  even  numbered  surfaces.  TRUNET  first  inter¬ 
polates  in  the  time  direction  bet%#een  corresponding  points  on 
the  two  surfaces  obtaining  the  data  at  points  on  a  constant 
time  plane.  It  then  equally  spaces  field  points  on  lines  con¬ 
necting  corresponding  shock  and  body  points  in  the  constant 
time  plane.  Flow  properties  at  these  new  field  points  are  ob¬ 
tained  by  linear  interpolation  between  three  adjacent  points 
in  the  constant  time  plane,  which  vere  previously  determined 
in  the  interpolation  in  the  time  direction.  Hence,  TRUNET 
calculates  two  new  surfaces  (an  even  nmbered,  and  an  odd 
numbered  surface),  but  also  interpolates  to  obtain  the  data 
on  a  constant  time  plane  at  equally  spaced  points. 
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TIMTER 


This  function  oMtrely  tpecifict  the  constant  time  at  which 
TRUMET  it  to  place  the  constant  time  plane.  This  time  should  be 
sosMwhere  bat%fean  the  two  even  numbered  surfaces  that  are  used 
in  the  interpolation  in  TRUNET.  The  time  is  usually  determined 
by  calculating  a  typical  point  on  the  new  surface  being  deter¬ 
mined.  If  the  time  specified  is  not  between  the  two  surfaces, 
extrapolation  rather  than  interpolation  will  be  carried  out  in 
TRUNET.  The  specification  of  the  time  was  removed  from  TRUMET 
to  the  small  subroutine,  TIMTER,  to  facilitate  easy  change  of 
the  specified  time. 


PATCH 

The  PATCH  subroutine  accepts  the  base  points  as  input, 
amd  determines  the  location  and  flow  properties  at  the  inter¬ 
mediate  base  points ,  so  that  the  domain  of  dependence  of  the 
difference  schMe  contains  the  domain  of  dependence  of  the 
partial  differential  equations.  This  insures  the  nimierical 
stability  of  calculation  (see  Appendix  B). 

PATCH  has  four  different  options  for  locating  the  inter¬ 
mediate  base  points  corresponding  to  values  of  1,  2,  3,  and  4 
of  its  first  arguMnt.  Its  second  argument  is  the  convergence 
test  value  used  in  the  iteration  loop.  The  first  option  locates 
the  intermediate  base  points  for  the  field  point  procedure,  as 
^■bown  in  Figure  20(a),  and  discussed  in  detail  in  Chapter  3. 

The  second  option  was  an  unsuccessful  scheme  for  the  body  and 
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•hock  point  proctdurea.  A  circle  was  fitted  ••  shown  in  Figure 
20(b),  end  the  points  of  tengency  between  the  field  points  end 
body  (or  shock)  points  were  used  es  the  location  of  the  Inter- 
■edlate  base  points.  This  scheme  did  ndt  completely  satisfy 
the  stability  criterion  and  a  mild  Instability  was  discovered 
In  the  results,  after  10  or  more  time  surfaces  had  been  cal¬ 
culated. 

The  third  and  fourth  options  are  successful  schemes  for 
the  shock  and  body  point  procedures,  respectively.  They  both 
use  the  pattern  shown  in  Figure  20(c).  First,  an  Intermediate 
shock  (or  body)  point  Is  determined  as  an  average  of  the  shock 
(or  body)  base  points.  The  Intermediate  base  points  are  then 
located  four  tenths  of  the  way  up  the  line  from  this  Inter¬ 
mediate  shock  (or  body)  point  to  the  field  base  points.  The 
only  difference  between  the  shock  and  body  point  options  Is 
that  the  intermediate  shock  point  Is  used  In  SHICPT,  while  the 
Intermediate  body  point  la  not  used  in  BOTPT. 

Linear  interpolation  Is  used  everywhere  in  PATCH.  This 
part  of  determining  the  location  of  the  Intermediate  base 
points  was  not  Included  In  the  FLDPT,  BDYPT  and  SHKPT  sub¬ 
routines,  so  that  If  higher  order  interpblatlon  were  desired 
at  a  later  tlsae.  It  could  be  easily  Included  by  merely  writing 
a  new  PATCH  subroutine. 

PATCH  has  three  error  checking  procedures  In  the  first 
two  options  where  iteration  Is  necessary.  If  the  Iteration  does 
not  converge  In  3O  cycles  through  the  loop,  the  calculation  Is 
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halted,  control  it  transferred  back  to  the  calling  program,  and 
an  error  statement  is  printed  off-line.  The  sets  of  linear 
algebraic  equations  are  solved  using  the  M.I.T.  Computation 
Center  library  routine  XSIMEQ.  If  over  or  underflow  occurs 
while  the  equations  are  being  solved  or  if  the  matrix  of  the 
equations  is  singular,  the  computation  is  halted  and  appro¬ 
priate  error  statements  are  printed  off-line.  XSIMEQ  is  used 
in  FLDPT,  BDTPT  and  SHKFT,  and  the  same  error  checking  procedure 
is  used  in  those  subroutines. 

FLDPT 

The  FLDPT  subroutine  calculates  the  coordinates  and  flow 
properties  of  a  new  point  in  the  field,  given  the  intermediate 
base  points  calculated  by  PATCH.  It  does  this  using  three 
iteration  loops;  two  for  the  third  order  Richmond  scheme  ,  to 
solve  the  Mach  cone  equations,  and  one  to  converge  the  solution 
to  the  compatibility  equations  to  second  order  in  the  step  sixe. 

The  convergence  tests,  as  proposed  in  Chapter  3i  have  the 

form 

>)_  .(n-i) 

^ _ Z _ ! _ <convergence  test  value 

suitable  normalizing  factor  (^>1) 

where  f  is  a  typical  quantity  to  be  converged,  and  the  in¬ 
equality  must  be  fulfilled  before  the  sequencing  is  released 
from  the  loop.  Typical  convergence  test  values  range  from 
0.0001  to  0.000001.  These  tests  are  employed  on  the  inner 
and  outer  loops,  but  it  was  found  that  the  middle  loop,  which 
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tcttf  the  coordinetes  of  the  new  point,  lOMtlMi  could  not 
converge  then  to  typicel  convergence  test  velues.  This  wes 
beceuse  the  reslduels  in  the  Mech  cone  equetions  had  already 
been  reduced  to  zero, to  eight  decimal  digits  which  is  the 
ultimate  accuracy  of  the  machine  using  single  precision 
arithnetic.  For  this  reason,  the  residuals  in  the  equations 
are  tested  in  the  middle  loop  tests  where  typical  convergence 
test  values  of  O.OOOOOO5  are  used. 

As  error  checks,  FLDPT  has  the  same  3O  cycle  limit  on  its 
looping  as  PATCH.  The  SHKPT  and  BDYPT  subroutines  also  have 
the  same  limits.  In  each  case,  control  is  returned  to  the  call¬ 
ing  program,  and  an  error  statement  is  printed  off-line  together 
with  the  z,  y,  and  t*  coordinates  of  the  point  which  did  not 
converge . 

SHKPT 

The  SHKPT  subroutine  calculates  the  coordinates  and  flow 
properties  of  a  new  point  on  the  shock  wave,  together  with  the 
components  of  the  shock  wave  normal  vector.  It  has  three  itera¬ 
tion  loops  as  in  the  FLDPT  subroutines.  The  error  checking  pro¬ 
cedures  and  convergence  tests  are  also  exactly  the  same  as  those 
in  FLDPT. 

It  was  found  that  in  using  equation  (3*30)  to  calculate 
ANst't  even  though  the  quantity  under  the  radical  did  not  con¬ 
verge  to  a  negative  number,  it  was  possible  for  it  to  take  on  a 
negative  value  during  the  outer  loop  iteration.  This,  of  course. 
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caused  the  aechlne  to  error  stop  beceuee  it  wee  not  progrenaied 
to  handle  the  iaaginary  nuaber.  This  trouble  was  eliminated 
by  testing  and  and  using  the  smaller  together  with 
in  the  expansion,  (3 >28).  Eliminating  the  largest  component 
and  solving  for  its  correction  using  equation  (3*30),  insures 
that  the  quantity  under  the  radical  will  not  turn  slightly 
negative  during  the  iteration.  This  procedure  has  been  in¬ 
corporated  in  the  programs  and  seems  to  eliminate  the  trouble. 


BDYPT 

This  subroutine  calculates  the  coordinates  and  flow 
properties  of  a  new  point  on  the  body  surface.  It  uses  the 
same  three  iteration  loops,  convergence  tests  and  error  checking 
procedures  as  the  SHKPT  and  FLDPT  subroutines  described  above. 

It  was  found  that  in  certain  expansion  regions  on  the  body 
surface ,  the  pressure  could  become  negative  in  the  outer  loop 
iteration  even  though  it  finally  converged  to  a  positive  quantity. 
This  would  cause  an  error  stop  in  the  ROEOET  or  ADET  functions 
where  the  pressure  is  raised  to  a  nonintegral  power.  To  eliminate 
this  trouble,  the  pressure  is  set  equal  to  zero  if  it  turns 
negative.  BDTPT  checks  the  final  value  of  the  pressure  before 
returning  control  to  the  calling  program,  and  prints  a  warning 
statement  off-line  if  the  final  converged  value  of  the  pressure 
is  zero. 
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All  r«f«r«nc«s  to  th«  function 

'dt  ^')  -  ^  (‘*•2) 

which  dosoribos  th«  body  aurfAce  have  been  removed  from  the 
BDTPT  aubroutine  and  placed  in  a  aubroutine  called  XCHAMG.  Thla 
waa  done  to  minimize  the  number  of  changea  neceaaary  to  conaider 
different  motiona  of  the  aame  body  or  in  certain  caaea ,  the  flow 
about  different  bodiea.  The  general  form  of  XCHARG  ia  given  in 
the  form  of  a  Hating  in  Appendix  C. 

XIKTER 

Thia  aubroutine  uaea  linear  Interpolation  to  determine  a 
flow  property  at  a  fourth  point  given  the  coordinatea  of  three 
other  pointa  at  which  the  flow  property  ia  known. 

ROEDET 

Thia  function  utilizea  the  perfect  gaa  equation 


where 

p  -  2117  Ib./ft.^ 
o 

-  0.002498  alug/ft.^ 
ar  - 1.4 

-  41.11 

R  -  0.06886  BTU/lbin.®R  -  1724  ft.^/aec^  °R 
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to  compute  the  density  given  the  pressure  and  entropy. 


ADET 

Given  Che  pressure  and  entropy,  this  function  determines 
the  speed  of  sound  using  the  perfect  gss  equation 


where  the  constants  are  given  above. 


THETA 

This  is  a  very  simple  function  i#hich  determines  an  angle 
in  radians  when  given  the  opposite  and  adjacent  sides  of  a  right 
triangle. 
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CHAPTER  5 

RESULTS  AND  CONCLUSIONS 


3.1  The  Feasibility  of  Calculating  Practical  Flow  Problems 
The  feasibility  of  carrying  out  a  practical  calculation 
of  an  unsteady  flow  field  hinges  on  whether  the  flow  can  be 
calculated  In  a  reasonable  number  of  time  steps.  An  estimate 
of  the  number  of  time  steps  required  can  be  obtained  in  the 
following  manner.  The  total  time  required  for  a  practical 

calculation  can  be  expressed  as 

^  I 

trrt  =  —  (5-1) 

4  ,  , 

Vave  average  particle  velocity  In  the  field,  L  Is  a 

typical  body  length  and  K  is  the  number  of  particles  which  m\ist 
pass  over  the  body  during  a  typical  transient  motion  or  dur¬ 
ing  one  cycle  of  an  oscillation.  K  can  also  be  considered 
the  number  of  body  lengths  a  particle  must  travel  during  a  ty¬ 
pical  transient  motion.  Now  the  average  time  step  Is  related 
to  the  space  mesh  size  by  the  following  well  known  estimate  of 


Its  order  of  magnitude. 


(5.2) 


This  can  be  easily  derived  by  considering  the  geometry  of  the 
characteristics.  Now  an  estimate  of  the  total  number  of  time 


steps  Is  given  by 


FKi 

-  / , 

1  j<_ 

— 

L 

(5.3) 


where 


£  - 


(5.^) 


Cl.gi^^  I  '  I® 

Now  ^can  probably  vary  between  the  extreme  cases  of 
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£  •  i 


M,  =  2(5 


L/ck  •  ZO  K  -  S 
which  gives  a  7"  of  10  to 

L  /ox  s  So  K  =  lO 
which  reaulta  In  a  oi  1000.  Hence,  It  appears  that  higher 
Mach  number  problems  which  do  not  require  too  fine  a  mesh 
spacing,  ^x/L,  can  be  attempted  with  this  method.  With  the 
present  programs,  which  utilize  linear  Interpolation,  the  upper 
limit  on  seems  to  be  approximately  50  because  of  growing 
errors  Introduced  by  the  Interpolation.  The  use  of  higher  order 
Interpolation  should  Increase  this  limit  by  a  factor  of  perhaps 
three  to  four. 

Another  area  where  the  solution  of  two-dimensional  un- 

hQ 

Steady  flow  arises  Is  In  hypersonic  small-dlsturbance  theory  . 
In  this  theory  the  flow  over  a  three-dimensional  body  Is  simpli¬ 
fied  to  consider  the  two-dimensional  unsteady  flow  In  plane 
slabs  which  sweep  over  the  body  at  the  free  stream  velocity. 

An  estimate  of  the  number  of  time  steps  required  to  do  this  sort 
of  problem  Is  given  by 

L 


7-  *  =  0 

At 


Jz.  A 


=  0 


(5.5) 


Even  for  l/Ax  »  50  and  «  5#  only  10  time  steps  are  necessary 
to  do  the  problem.  Hence,  It  appears  that  most  problems  utiliz¬ 
ing  hypersonic  small-dlsturbance  theory  can  be  easily  handled. 


5.2  Operational  Results  for  the  Programs 

Hie  results  of  a  typical  field  polr.c  calculation  are  pre¬ 
sented  In  Table  I.  These  results  from  PLDPT  subroutine 
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are  typical  In  the  amount  of  looping  required,  the  rate  of  con¬ 
vergence  and  the  accuracy  obtained  In  the  Iterations.  Similar 
results  are  obtained  with  the  BDYPT  and  SHKPT  subroutines.  It 
should  be  noted  that  the  Inner  loop  usually  requires  three 
cycles  of  the  Iteration  while  the  middle  loop  usually  satis¬ 
fies  the  convergence  criterion  the  first  time  through  the  loop. 
I^e  outer  loop,  which  In  general  governs  the  amount  of  computa¬ 
tion  time  required,  usually  takes  four  to  six  cycles  of  the 
Iteration  to  converge.  The  results  are  given  to  eight  digits 
In  Table  I  because  this  Is  the  number  of  digits  used  by  the 
machine  In  the  calculation.  Only  the  first  four  or  five  digits 
are  significant,  however. 

It  can  also  be  seen  from  Table  I  that  the  calculation  Is 
converging  very  rapidly  for  at  the  bottom  of  the  table  are  pre¬ 
sented  the  results  of  the  same  calculation  with  smaller  conver¬ 
gence  test  values.  These  results  Indicate  that  the  Inner  loop, 
for  example,  has  actually  been  converged  to  a  convergence  test 
value  one  order  of  magnitude  smaller  than  required  In  the  first 
calculation.  An  additional  two  cycles  through  the  outer  loop 
satisfies  a  convergence  test  value  two  orders  of  magnitude  small¬ 
er  than  the  original  value. 

These  results  might  Indicate  that  the  convergence  test 
criteria  could  be  eliminated  and  the  calculation  carried  out  with 
a  fixed  set  of  cycles  through  each  loop.  This,  however.  Is  not 
recommended  without  a  very  thorough  study  of  the  ranges  of 
values  of  all  Ir^ut  parameters.  It  Is  very  possible  that  for 
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certain  areas  of  flow  fields  or  certain  network  spaclngs  addi¬ 
tional  Iteration  might  be  necessary  to  obtain  reasonable  accur¬ 
acy. 

The  computation  time  required  to  determine  the  coordinates 
and  flow  properties  at  a  new  point  Is  very  difficult  to  deter¬ 
mine  precisely.  It  depends  on  the  convergence  test  values  and 
on  the  variation  of  the  flow  properties  at  the  base  points  (l.e. 
the  size  of  the  net  spacing) .  Estimates  of  the  time  required 
to  compute  one  point  have  been  obtained  In  the  following  man¬ 
ner.  The  total  time  required  to  read  the  data  from  the  tape, 
calculate  two  new  time  surfaces  with  a  total  of  600  to  1000 
points  and  vrrlte  the  data  back  on  tape  was  divided  by  the  total 
n\imber  of  points  calculated.  In  general  the  time  ranged  from 
0.25  to  0.45  seconds  per  point  with  an  average  approximately 
In  the  middle.  This  time,  of  course.  Includes  body  and  shock 
points  and  also  the  time  required  to  load  the  Input  data  Into 
the  proper  subroutines.  It  Is,  however,  a  conservative  estimate 
of  the  time  required  to  determine  a  new  field  point. 

The  executive  program  and  subroutines  listed  In  Appendix 
C  together  with  required  Input-output,  library,  and  error 
tracing  routines  take  up  about  20,000  words  of  the  magnetic 
core  storage.  This  does  not  Include  the  common  storage  area 
where  all  the  data  Is  stored.  This  storage  requirement  can 
be  reduced  by  Improving  the  subroutines  as  discussed  In  Sec¬ 
tion  6.1.  Several  error  checking  and  debugging  routines  brought 
In  from  the  library  could  also  be  eliminated.  It  Is  felt  that 
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the  minimum  storage  required  for  the  programs  lies  somewhere 
between  10,000  and  15,000  words,  perhaps  closer  to  the  larger 
figure.  Thus,  the  original  estimate  In  Chapter  4  that  about 
half  of  the  32,000  word  core  storage  would  be  needed  for  the 
storage  of  programs  seems  correct. 

The  PLDPT,  BDYPT,  and  SHKPT  subroutines  were  Individually 
checked  Initially  using  Input  data  obtained  from  analytical 
solutions  to  rather  simple  flow  fields.  The  output  was  then 
checked  against  the  known  analytical  solutions.  The  one-dlmen- 
slonal  unsteady  motion  of  a  piston  In  a  cylinder  of  gas  was 
used  to  check  the  PLDPT  and  BDYPT  subroutines.  The  motion  of 
normal  and  oblique  plane  shock  waves  was  used  to  check  SHKPT. 
These  checks  were  not  exhaustive  In  that  they  covered  large 
ranges  of  the  Input  parameters.  Due  to  time  limitations  the 
checks  were  carried  out  only  until  reasonable  assurance  was 
obtained  that  the  subroutines  were  working  properly.  The  sub¬ 
routines  were  given  additional  checks  In  calculating  complete 
flow  fields  and  the  results  of  these  checks  are  given  In  Sec¬ 
tion  5.4. 

5.3  Difficulties  Encountered  and  Their  Solutions 

In  the  course  of  programming  and  debugging  this  method  of 
characteristics  several  difficulties  were  encountered  which  had 
to  be  overcome  or  avoided  In  order  to  carry  out  practical  calcu 
latlons.  It  Is  felt  that  a  brief  description  of  these  dlfflcul 
ties  might  help  those  who  attempt  to  utilize  or  extend  the  pro¬ 
cedures  presented  here. 


93 


Probably  the  moat  fundamental  difficulty  encountered  was 
the  numerical  Inatablllty  which  aroae  when  the  tetrahedral  char¬ 
acteristic  line  network  waa  used  Initially.  This  trouble  was 
eliminated  by  gaining  an  understanding  of  the  basic  cause  of 
the  Instability  (see  Appendix  B)  and  then  utilizing  this  under¬ 
standing  to  propose  a  modification  to  the  network  such  that 
stability  was  obtained.  The  modified  tetrahedial  characteristic 
line  network  was  the  solution  to  this  difficulty  and  the  PATCH 
subroutine  was  the  programmed  form  of  the  solution. 

It  was  pointed  out  In  Chapter  4  that  the  location  of  a  new 
point  Is  a  function  of  the  spacing  of  the  base  points  In  the 
Initial  surface,  the  local  flow  velocity  and  the  local  speed  of 
sound.  The  last  two  factors  cannot  be  controlled  so  that.  If 
no  effort  Is  made  to  control  the  positions  of  the  base  points, 
the  net  cam  become  skewed.  This  In  fact  Is  what  happened  In  all 
the  flow  fields  calculated.  After  several  time  surfaces  had  been 
determined  the  last  surface  on  which  the  flow  was  know  became 
skewed  so  that  It  was  no  longer  close  to  being  a  constamt  time 
plane.  When  the  skewing  got  very  bad  the  approximations  which 
are  used  to  determine  the  first  estimate  of  the  coordinates  of 
the  new  point  give  very  poor  estimates.  The  Iteration  process 
to  determine  the  coordinates  of  the  points  then  takes  a  large 
number  of  cycles  to  converge  and  eventually  Is  unable  to  converge 
at  all.  Hence,  In  order  to  be  able  to  continue  the  calculation, 
the  skewing  must  be  controlled  or  eliminated.  The  only  way  to 
do  this  Is  to  control  the  location  of  the  points  In  the  Initial 
surface.  This  was  done  by  using  Interpolation  periodically  to 
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eliminate  the  net  skewing.  T^e  TRUNET  and  ADDROW  subroutines 
described  In  Chapter  U  were  written  to  perform  this  interpolation 
and  to  redistribute  the  points  In  the  Initial  surface. 

Some  difficulty  was  encountered  In  choosing  the  proper 
form  of  the  convergence  tests  In  the  three  loops.  As  discussed 
In  Chapter  h.  It  was  found  that  for  some  of  the  smaller  net  spac- 
Ings  the  residuals  In  the  Mach  cone  equations  were  reduced  to 
zero  to  eight  decimal  digits  before  the  coordinates  had  converged 
to  four  or  five  digits.  Hence,  It  was  found  necessary  to  utilize 
residual  tests  In  the  middle  loops  of  the  PLDPT,  SHKPT,  and  BDYPT 
subroutines.  The  outer  loop  tests  also  caused  some  trouble 
In  these  subroutines  when  one  of  the  velocity  components  was 
close  to  zero.  This  trouble  was  eliminated  by  testing  either 
the  pressure  or  the  total  velocity  magnitude,  whichever  was 
greater.  In  this  way  there  Is  almost  no  chance  of  trying  to 
test  a  quantity  whose  magnitude  lo  almost  zero.  Quantities 
which  pass  through  zero  are  difficult  to  use  In  convergence 
tests  because  they  cannot  be  suitably  normalized. 

3.^  The  Results  for  Complete  Flow  Fields 
Blunt  Body  Flows 

The  flow  over  a  circular  cylinder  with  Its  axis  normal  to 
a  supersonic  free  stream  has  been  considered  as  an  Interesting 
and  challenging  complete  flow  field  to  study.  'Rils  flow  Is  of 
Interest  In  that  It  has  a  detached  shock  wave  which  bounds  the 
field  on  one  side  and  the  body  surface  which  forms  the  boundary 
on  the  other.  The  flow  field  Is  challenging  because  It  contains 
both  subsonic  and  supersonic  regions  and  hence  also  has  large 
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variations  In  the  flow  properties.  Obtaining  the  Initial  steady 
flow,  from  which  unsteady  flows  can  be  calculated.  Is  by  no  means 
a  simple  matter.  The  steady  flow  Involves  the  solution  of  a 
mixed  elliptic-hyperbolic  problem  with  a  free  boundary.  Belot- 

liq 

serkovskll  ^  has  solved  this  steady  flow  numerically  using  the 
method  of  Integral  relations. 

Figure  21  shows  the  Initial  data  net  used  for  the  blunt 
body  flows  presented  here.  "Hiese  points  and  the  flow  properties 
at  each  of  the  points  were  taken  from  the  tabular  results  pre¬ 
sented  by  Belotserkovskll .  "nie  net  was  made  more  uniform  by  add¬ 
ing  four  rings  of  points  by  linear  Interpolation.  The  properties 
were  given  to  three  or  four  digits  In  the  tabular  presentation. 

Two  strips  were  used  In  the  method  of  Integral  relations  to  ob¬ 
tain  this  initial  flow.  A  complete  symmetric  net  was  obtained 
by  reflecting  the  points  about  the  x-axls. 

The  highest  free  stream  Mach  number  for  which  Belotser¬ 
kovskll  presented  results  was  five.  Unfortunately,  when  (5.^) 

Is  applied  to  this  low  Mach  number  case  (with  K  -  3^  £  -  0.1 
and  ^x/L  >  0.03)  the  estimated  total  number  of  time  steps  required 
for  a  typical  problem  Is  200.  Obviously  this  Is  not  a  hyper¬ 
sonic  flow  and  a  complete  calculation  would  require  an  excessive 
amount  of  conputatlon  time.  For  example,  with  the  net  of  330 
points  considered  here,  200  time  steps  at  one  third  of  a  second 
per  point  would  require  six  and  one  half  hours  of  computer  time. 
Even  If  this  large  amount  of  computer  time  were  available,  the 
linear  Interpolation  utilized  In  these  Initial  programs  would 
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Introduce  an  Increaalng  Inaccuracy  which  would  make  the  solution 
obtained  after  as  many  as  200  time  steps  only  qualitatively  cor¬ 
rect.  However,  because  of  the  unavailability  of  hilgher  Mach 
number  solutions  of  two-dimensional  perfect  gas  blunt  body  flows, 
the  M^  =  5  data  of  Belotserkovskll  can  be  used  to  check  the 
feasibility  of  carrying  out  practical  calculations  by  consider¬ 
ing  at  least  the  Initial  part  of  some  unsteady  motions.  TTie 
results  obtained  starting  with  this  Initial  data  and  considering 
several  types  of  unsteady  body  motion  are  discussed  In  the  fol¬ 
lowing. 

Pour  different  blunt  body  flow  cases  were  calculated  for 
presentation  here.  A  steady  flow  was  calculated  by  holding  the 
body  surface  unchanged  with  time.  The  second  case  warped  the 
body  symmetrically  from  a  circular  cross-section  to  an  elliptic 
cross-section  In  a  short  length  of  time.  A  similar  case  was  nan 
with  the  warping  done  asymmetrically  so  that  the  final  elliptic 
section  was  at  an  angle  of  attack  of  20°.  "nie  fourth  case  oscil¬ 
lated  the  body  vertically  (In  the  y-dlrectlon)  nflth  a  1  -  cos  oat 
time  dependence. 

* 

The  steady  flow  was  run  to  ascertain  how  well  the  calculation 
procedure  maintained  the  Initial  steady  flow.  It  vras  expected 
that  the  linear  Interpolations  in  PATCH,  TRUNET  and  ADDROW 
would  Introduce  growing  systematic  errors.  As  might  be  ejqjected, 
the  flow  field  tended  to  be  "smoothed”  as  the  calculation  pro¬ 
ceeded  In  time.  When  linear  Interpolation  Is  used  to  repeatedly 
approximate  a  function,  the  points  which  describe  the  function 
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tend  to  move  toward  the  center  of  curvature  of  the  function  as 
■hown  In  the  following  simplified  sketch. 


o  Initial  points 

X  after  first  linear 
Interpolation 
•  after  second  linear 
Interpolation 

curvature 

The  rate  at  which  they  move  is  proportional  to  the  curvature. 

Of  course,  the  numerical  solution  of  the  partial  differential 
equations  considered  here  Is  considerably  more  complicated,  but 
the  linear  interpolation  might  be  expected  to  have  the  same 
qualitative  effect.  This  can  be  seen  In  Plg\u?e  22  where  the 
body  surface  pressure  Is  plotted.  Note  that  the  origin  has  been 
suppressed  on  the  ordinate  of  the  plot. 

The  time  coordinates  of  the  time  surfaces  for  the  four 
blunt  body  cases  are  given  In  Table  II.  When  only  one  time 
Is  given  TRUNET  was  used,  and  when  a  range  of  time  Is  given  the 
regular  executive  program  was  used. 

The  body  surface  pressure  was  found  to  be  the  quantity 
which  changed  the  most  rapidly  with  time.  The  velocity  conq;>onents 
as  shown  In  Figures  23  and  24  did  not  change  as  rapidly  In  time. 

In  fact  they  drifted  only  slightly  from  their  Initial  values. 

The  pressure  on  the  middle  ring  of  field  points  shown  In  Figure 
23  has  the  same  smoothing  tendency.  Note  that  the  ciu*vature  In 
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Figure  25  Is  less  than  the  curvature  In  Figure  22  on  the  body 
surface  and  that  the  rate  at  which  the  pressure  falls  with  time 
Is  less  than  that  In  Figure  22.  Also  note  that  In  all  flow 
properties  except  the  u  component  of  velocity,  which  changes 
only  slightly,  the  tendency  Is  for  the  flow  property  to  move 
toward  Its  center  curvature  as  would  be  expected  If  the  linear 
Interpolation  were  causing  the  drift. 

In  Figure  26  the  pressure  along  the  axis  of  symmetry  from 
the  body  to  the  shock  Is  plotted  versus  the  x-coordlnate. 

Again  note  that  the  origin  on  the  ordinate  has  been  extremely 
suppressed  on  this  plot.  It  appears  that  the  falling  pressure 
on  the  body  surface  caused  by  the  linear  Interpolation  and  the 
large  "curvature"  In  pressure  on  the  body  surface  actually  In¬ 
troduces  a  mild  expansion  wave  Into  the  field.  The  drift  due 
to  the  Interpolation  appears  to  reach  such  a  magnitude  that  It 
Is  "picked  up"  suid  propagated  through  the  field  as  an  expansion 
wave. 

The  data  from  the  two  strip  method  of  Integral  relations 
used  here  as  Initial  data  should  be  accurate  to  at  least  three 
digits.  The  two  strip  results  are  quite  good  as  compared  to  the 
one  strip  results  frequently  used,  which  can  be  as  much  as  10^ 
or  more  too  low  for  the  pressure  on  the  body  surface  In  the  expan 
Sion  region  beyond  the  sonic  llne^^.  Thus,  It  cannot  be  reasoned 
that  the  unsteady  method  of  characteristics  Is  trying  to  adjust 
the  Initial  steady  flow  data  to  a  more  accurate  solution. 

Perhaps  a  smaller  net  spacing  would  have  reduced  the  drift 
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due  to  the  linear  Interpolation  (a  abort  teat  caae  Indicated  that 
It  did),  but  becauae  the  finer  net  reaulta  In  a  shorter  time 
atep  (aee  equation  (5.3)  )*  the  accuracy  In  calculating  a  given 
length  of  time  for  the  flow  would  probably  be  about  the  same. 
Hence,  It  appeara  that  higher  order  Interpolation  will  be  required 
to  Improve  the  accuracy  of  the  calculation.  The  particular  net 
size  choaen  waa  alao  a  compromlae  to  keep  the  total  computer 
time  required  to  carry  out  50  time  atepa  for  each  of  the  four 
blunt  body  caaea  at  a  reasonable  value.  Each  case  required  about 
one  and  one  half  hours  of  IBM  709^  computation  time  to  calcu¬ 
late  the  50  steps. 

It  can  be  seen  In  Figures  22  and  25  that  a  small  amount 
of  asymmetry  has  developed  In  the  field  at  the  latest  time  steps. 
This  was  caused  by  a  minor  asymmetry  In  the  organization  and 
sequencing  of  the  calculation  by  the  executive  program  and  TRUNET 
and  can  probably  be  eliminated  by  slight  modifications  to  these 
programs.  It  was  planned  to  call  TRUNET  at  alternating  time 
steps  and  to  call  ADDROW  at  every  fourth  step,  but  as  can  be 
seen  from  Table  II  this  sequencing  was  not  followed  exactly. 

The  sequence  was  disrupted  because  of  human  operational  errors 

\ 

In  the  running  of  the  programs.  Because  of  the  large  amount  of 
conqputer  time  Involved,  the  cases  were  not  rerun  to  follow  the 
sequencing  exactly  and  some  asymmetry  In  the  final  data  resulted. 
It  will  also  be  noted  that  the  range  of  y  over  which  the  calcu¬ 
lation  was  carried  out  decreased  slightly  as  the  calculation 
was  continued  In  time.  This  was  because  the  Initial  data  did 
not  continue  out  Into  the  supersonic  region  far  enough  so  that 
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the  Mach  cones  leaned  far  enough  downstream  to  make  their  Inter¬ 
sections  lie  even  further  downstream.  TTius,  the  total  area  cover¬ 
ed  by  the  calculation  on  each  time  surface  did  not  Increase,  but 
decreased  at  such  a  slow  rate  that  It  did  not  Introduce  much 
of  a  restriction  In  the  four  blunt  body  cases  considered. 


A  typical  fluid  particle  entering  the  shock  wave  at 
about  30°  above  the  axis  of  symmetry  (near  y  »  0.0  In  Figure  21) 
has  an  average  velocity  through  the  field  of  2000  fps.  The  total 
time  calculated  In  the  50  time  steps  was  about  0.00025  seconds. 
The  particle  moved  0.5  feet  In  the  50  time  step  and  passed 
about  half  way  through  the  flow  field  being  calculated.  In 
order  to  watch  five  such  particles  pass  through  the  field  500 
time  steps  would  have  to  be  calculated.  This  agrees  with  the 
order  of  magnitude  estimate  made  above  for  this  particular  net. 

The  second  2Lnd  third  blunt  body  cases  considered  were 


specified  by  the  following  equations  for  the  body  surface: 

for  0  <  t'  4  UB  £ 

-h  2p(  '1^0  (5.6a) 


and  for  t'>  U  B 


(&()(,  =  X  j  f  <  (t  -  {.j  3ifi  2  0^ 


Where 


**■«  -!  -  o 


(5.6b) 


/{  -  I.OS 

yp  -  o.ife  tt 

andO(*  0°  for  the  symmetric  case  andO(=  20®  for  the  asymmetric 


case  which  results  In  an  ellipse  whose  semi-major  axis  Is  55^ 


longer  than  Its  semi-minor  axis.  The  semi-major  axis  lies  along 
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the  x-axl8  for  the  symmetric  case  and  has  a  positive  angle  of 
attack  of  20°  for  the  asymmetric  case.  Note  that  the  transition 
from  the  circle  to  the  ellipse  takes  place  8nx>othly  using  one 
half  period  of  a  cosine  function. 

Figure  27  shows  the  body  surface  pressure  for  the  sym¬ 
metric  wazplng  case.  The  pressure  rises  to  a  maximum  on  the 
eighth  time  surface  when  the  body  surface  has  Its  maximum  velo¬ 
city.  It  then  falls  to  a  level  slightly  below  the  Initial  pres¬ 
sure  due  to  the  slowing  of  the  body  motion.  On  the  sixteenth 
surface  the  body  motion  has  stopped.  The  pressure  then  tends  to 
drift  doimward  as  In  the  steady  flow  In  Figure  22.  It  appears 
that  the  steady  flow  downward  drift  has  been  superimposed  on  the 
basic  unsteady  flow  and  that.  If  the  drift  were  removed  with 
higher  order  Interpolation,  the  stagnation  point  pressure  would 
return  to  the  same  value  It  had  Initially  as  Is  to  be  expected. 

The  compression  wave  traveling  from  the  body  surface 
toward  the  shock  wave  can  be  seen  In  Figure  28.  The  peeik  pres¬ 
sure  In  '  e  wave  falls  as  the  wave  moves  away  from  the  body 
because  of  the  two-dimensional  nature  of  the  wave.  The  wave 
speed  relative  to  the  fixed  coordinate  system  as  determined  from 
the  plot  Is  approximately  2500  fps.  Initially  and  about  I500  fps. 
at  the  later  time  surfaces.  Due  to  the  low  fluid  velocity  In  the 
stagnation  region,  these  values  are  approximately  equal  to  the 
speed  of  sound  In  the  field  at  these  locations.  For  example, 
the  speed  of  sound  near  the  body  surface  Is  approximately  2200 
fps.  so  that  a  compression  wave  only  slightly  stronger  than  a 
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Bound  wave  has  been  obtained.  Note  that  a  small  expansion  region 
has  developed  In  front  of  the  compression  wave.  Of  course  this 
could  not  have  been  Introduced  by  the  body  motion,  but  appears 
to  arise  from  the  pressure  decrease  caused  by  linear  Interpola¬ 
tion  which  was  shown  In  Figure  26  In  the  steady  flow  eind  dis¬ 
cussed  above. 

It  can  be  seen  from  Figure  28  that  the  net  Is  very  coarse 
In  the  x-dlrectlon.  It  Is  almost  too  coarse  to  Indicate  pro¬ 
perly  the  structure  of  the  wave.  It  was  necessary,  however,  to 
use  this  coarse  a  net  and  a  fast  body  motion  In  order  to  carry 
out  the  calculation  In  a  reasonable  length  of  computer  time. 

In  the  48  time  steps  calculated  the  wave  has  Just  reached 
the  shock  wave.  The  shock  wave  reacts  properly  to  the  wave 
In  that  It  moves  slightly  toward  the  body  when  the  small  expan¬ 
sion  region  first  hits  and  then  moves  more  rapidly  away  from 
the  body  when  the  stronger  compression  wave  strikes.  From  the 
decrease  In  amplitude  of  the  wave  In  going  from  the  body  to  the 
shock  once.  It  appears  that  the  transient  time  required  for  the 
flow  to  attain  a  steady  state  might  be  four  times  what  has  been 
calculated  here.  This  agrees  with  the  crude  estimate  of  200 
time  steps  given  at  the  first  of  this  section. 

The  results  of  the  asymmetric  warping  case  are  shown  In 
Figures  29  and  30.  The  body  surface  pressure  rises  and  falls 
In  exactly  the  same  way  as  in  the  symmetric  case,  but  the  dis¬ 
tribution  of  the  pressure  along  the  body  Is  now  asymmetric.  On 
the  eighth  time  surface  where  the  peak  pressure  Is  attained  the 


103 


stagnation  point  has  shifted  up  to  a  point  9°  above  the  original 
axis  of  symmetry.  After  the  body  warping  has  stopped  the  stag¬ 
nation  point  holds  fairly  steady  at  5°. 

« 

The  velocity  vectors  on  a  portion  of  the  l6th  and  3^th 
time  surfaces  jire  shown  In  Figure  30.  An  Idea  of  how  far  the 
compression  wave  has  moved  In  the  field  by  the  l6th  and  3^th  time 
surfaces  can  be  obtained  from  Figure  28.  In  the  l6th  time  sur¬ 
face  of  Figure  30  (a)  the  compression  wave  propagating  toward  the 
shock  from  the  body  can  be  clearly  seen.  Nc te  that  a  region  In 
which  the  gas  Is  actually  flowing  upstream  has  developed  next 
lo  the  body.  The  asymmetry  Is  evident  In  the  compression  wave 
between  the  wave  front  and  the  body.  This  asymmetry  can  also 
be  seen  In  Figure  30  (b)  where  the  compression  wave  has  traveled 
three  quarters  of  the  distance  from  the  body  to  the  shock.  The 
compression  wave  has  spread  and  weakened  so  much  that  there  Is 
no  longer  an  upstream  flow  region.  The  approximate  location  of 
the  streamline  which  passes  through  the  stagnation  point  has  been 
sketched  on  the  figure. 

For  the  fourth  blunt  body  case,  the  vertical  oscillation 
Is  described  by  a  body  surface  equation  of  the  form 

AjcosBt'- ijY’-i  -  O  (5.7) 

where 

A  -  0.02  ft. 

B  -  25  ft."^ 

The  frequency  of  the  motion  Is  7960  cycles  per  second  and  the 
reduced  frequency  Is 

k  =  W^/Uo  =  I0.2&  (5.8) 
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Thl8  high  frequency  was  chosen  so  that  several  periods  of  mo¬ 
tion  could  be  calculated  In  the  5^  time  steps  considered. 

The  results  of  the  calculation  In  the  form  of  body  sur¬ 
face  pressure  plots  are  given  In  Figure  31.  The  quasi-steady 
prediction  of  the  maximum  amplitude  pressure  distribution  Is 
shown  In  Figure  31(a).  This  quasi-steady  distribution  merely 
shifts  Instantaneously  adjusting  to  the  new  effective  free 
stream  velocity  obtained  by  vectorlally  subtracting  the  body 
velocity  vector  from  the  steady  free  stream  velocity  vector.  It 
can  be  seen  In  the  figure  that  the  stagnation  point  pressure  on 
the  fourth  and  sixth  time  surfaces,  which  fall  on  either  side 
of  the  maximum  amplitude  time,  agree  rather  well  with  the  quasi- 
steady  value.  The  location  of  the  stagnation  point  Is  at  about 
a  50%  greater  y-coordlnate  than  predicted  by  the  quasi-steady 
result.  Itiere  Is  also  a  marked  asymmetry  In  the  pressure 
distributions  which  Is  not  predicted  by  the  quasi-steady  re¬ 
sults. 

No  shift  In  the  phase  angle  between  the  pressure  oscil¬ 
lation  and  the  motion  of  the  body  could  be  detected  in  the  two 
cycles  of  motion  calculated.  The  stagnation  point  pressure  at 
the  maximum  pressure  amplitude  positions  fell  as  the  calculation 
was  continued  In  time.  This  Is  probably  the  same  effect  noted 
In  Figure  22  for  the  steady  flow  because  the  stagnation  point 
pressure  Is  falling  at  exactly  the  same  rate  of  about  IW 
per  time  step.  The  Introduction  of  higher  order  Interpolation 
should  eliminate  this  downward  drift. 

Wedge  Flow  Check  Cases 

The  first  complete  flow  field  to  be  euialyzed  was  super¬ 
sonic  flow  over  a  wedge  of  a  perfect  gas  with  Y  •  1.^.,  Itie 
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executive  program  wae  a  simplified  form  of  the  executive  used 
for  the  blunt  body  and  described  In  Chapter  4.  Net  conflgiira- 
tlons  similar  to  the  one  shown  In  Figure  32  were  used  In  approx¬ 
imately  230  separate  flow  field  calculations. 

The  main  purpose  of  the  calculations  was  to  check  the 
various  parts  of  the  programs  against  known  solutions.  No 
attempt  was  made  to  carry  out  useful  flow  calculations.  Note 
In  Figure  32  that  no  attempt  was  made  to  extend  the  calculation 
to  the  point  of  contact  of  the  shock  wave  with  the  wedge  tip. 

'nils  would  require  a  separate  program  to  calculate  the  special 
point  which  lies  both  on  the  shock  wave  and  the  wedge  surface, 
'nils  attached  shock  point  could  bj  easily  calculated  on  a  quasi- 
steady  basis  anyway.  If  It  were  needed.  Because  the  flow  Is 
entirely  supersonic,  the  area  which  can  be  calculated  from  the 
initial  data  decreases  with  each  time  step.  One  column  of  points 
running  from  the  shock  to  the  wedge  is  lost  with  each  time  step. 
Hence,  only  as  many  time  surfaces  as  there  are  columns  of  points 
can  be  calculated.  This  limits  the  extent  of  the  calculation,  but 
this  simple  test  case  was  only  meant  to  check  the  complete  pro¬ 
cedure. 

This  flow  field  served  as  a  very  good  check  case.  Various 
programming:  errors  became  apparent  and  mai^^rogram  improvements 
were  determined.  One  major  discovery  was  the  numerical  Insta- 

bility  discussed  above  in  Section  3*3*  This  instability  was 

J 

studied  and  eliminated  using  this  wedge  flow. 

A  few  Mach  three  flows  were  calculated,  but  most  of  the 
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flowB  had  a  free  stream  Mach  number  of  nine  and  a  30°  shock 
wave  angle.  Hundreds  of  runs  were  made  either  with  the  wedge 
surface  held  steady  In  time  or  Impulsively  moving  It  up  or  down 
at  right  angles  to  the  free  stream  velocity.  A  few  runs  were 
made  with  changing  free  stream  conditions  or  with  an  Irregularly 
shaped  shock  wave.  Most  of  these  calculations  were  made  while 
the  programs  were  being  debugged  and  Improved  so  that  only  a  few 
of  the  last  calculations  performed  can  be  considered  meaningful 
results . 

It  was  felt  that  rather  than  repeat  many  of  the  wedge  flow 
calculations  after  the  programs  worked  well,  more  could  be  learned 
by  proceeding  to  more  Interesting  and  complicated  flow  fields. 
Hence,  only  the  results  of  the  last  few  wedge  flow  calculations 
are  presented  here.  When  the  wedge  surface  was  held  steady,  the 
steady  flow  results  held  constant  to  six  and  seven  digits.  This 
was  to  be  expected  because  the  flow  properties  were  constant  In 
the  field  between  the  shock  and  the  wedge  surface.  No  outer  loop 
Iteration  was  necessary  and  only  the  coordinates  of  new  points 
had  to  be  determined.  The  results  of  four  calculations  In  which 
the  wedge  surface  was  Impulsively  set  Into”  motion  upward  at 
t^  =  0  are  presented  In  Figure  33.  The  wedge  was  moved  upward 
at  four  different  velocities.  This  gave  rise  to  a  compression 
wave  which  propagated  from  the  wedge  s\irface  out  toward  the  shock 
wave.  This  Impulsive  upward  velocity  Is  equivalent  to  an  In¬ 
stantaneous  Increase  In ^the  wedge  angle.  The  Instantaneous 
acceleration  also  corresponds  to  a  comer  or  discontinuity 
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on  the  body  surface.  The  flow  about  bodies  with  discontinuous 
surfaces  or  comers  cannot  be  calculated  with  the  present  pro¬ 
grams  as  was  pointed  out  In  Chapter  3*  Small  comers  should 
Introduce  only  small  errors  In  the  flow  calculation,  but  these 
errors  will  grow  with  the  size  of  the  discontinuity.  The  results 
of  these  wedge  flows  can  be  used  to  determine  the  size  of  the 
body  surface  discontinuity  which  can  be  tolerated  without  Intro¬ 
ducing  large  errors  Into  the  flow  calculation. 

In  Figure  33(a)  the  pressure  on  the  wedge  surface  for 
steady  flow  calculated  with  the  oblique  shock  equations  (Refer¬ 
ence  No.  50)  Is  plotted  versus  the  pressure  on  the  wedge  surface 
at  the  latest  time  that  could  be  calculated  before  no  more  points 
were  available  to  continue  the  calculation.  The  disagreement  of 
the  pressures  Increases  with  the  Instantaneous  wedge  angle  change 
for  two  reasons.  One  Is  that,  for  larger  Impulsive  angle  changes 
which  correspond  to  body  surface  discontinuities,  the  finite 
difference  approximations  become  less  accurate  as  mentioned 
above.  Actually,  very  large  angle  changes  could  lead  to  the 
formation  of  shock  waves.  A  tendency  for  the  compression  waves 
to  steepen  was  noted  In  some  of  the  calculations  even  though 
the  formation  of  a  shock  cannot  be  accurately  calculated  with 
the  present  programs.  The  second  reason  for  the  Increasing  dis¬ 
agreement  Is  that  the  flow  field  had  not  settled  to  a  completely 
steady  state  at  the  time  of  the  last  step  In  the  calculation. 

It  appears  that  the  body  surface  discontinuity  Is  too  large  for 
the  fourth  run  which  has  an  Instantaneous  angle  change  of  2.37°. 
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Hence,  the  area  to  the  right  of  the  shaded  region  probably  corres¬ 
ponds  to  a  body  surface  discontinuity  which  Is  too  large. 

Figure  33(b)  compares  the  Mach  number  behind  the  shock 
from  the  oblique  shock  equations  vrlth  the  Mach  number  from  the 
numerical  calculations.  The  Mach  number  was  obtained  Indirectly 
from  the  numerical  results  by  using  p  and  s  to  determine  a  which 
was  then  divided  Into  the  velocity  magnitude.  These  Mach  number 
results  seem  to  be  a  little  less  accurate  than  the  pressure  re¬ 
sults,  but  seem  to  have  the  same  trend. 

Figure  3^  shows  a  typical  time  history  of  the  pressure 
along  the  wedge  surface.  These  results  are  from  case  number  3 
of  Figure  3^*  The  variation  of  pressure  along  the  surface  on  the 
8th,  12th,  and  l8th  time  surface  Is  well  within  the  accuracy 
to  be  expected  for  the  convergence  test  values  of  0.001  used  In 
the  calculation.  The  variation  on  the  4th  surface  Is  slightly 
larger  than  might  be  expected,  but  It  Is  possible  that  the  abrupt 
nature  of  the  Instantaneous  vertical  motion  caused  some  Initial 
spurious  transient.  It  Is  also  possible,  however,  that  the  vari¬ 
ation  results  from  a  combination  of  the  organization  used  In 
the  executive  program  and  the  Iterations  In  the  various  subrou¬ 
tines  so  that  the  results  are  not  as  accurate  as  the  single  con¬ 
vergence  test  value  of  0.001  might  Indicate. 

3 . 3  Conclusions 

The  solution  of  unsteady  two-dimensional  flow  fields  by 
the  method  of  characteristics  Is  feasible  and  practical  at* 
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the  present  time  utilizing  presently  available  high  speed  elec¬ 
tronic  digital  computers,  "nie  solution  of  the  more  difficult 
problems  Involving  nonequilibrium  thermodynamics  or  three-dimen¬ 
sional  unsteady  flows  Is  marginal  with  existing  computers  but 
should  become  completely  feasible  with  the  next  generation  of 
computers  which  will  appear  shortly. 

The  problems  being  considered  here  Involve  large  amounts 
of  Input  and  output  data.  For  example,  the  output  data  from  the 
four  blunt  body  cases  discussed  In  the  previous  section  was  In 
the  form  of  approximately  one  half  million^ eight  digit  numbers. 

The  handling  and  analysis  of  this  data  was  by  no  means  a  simple 
task.  The  development  of  methods  and  devices  for  quickly  hand¬ 
ling  large  amounts  of  data  and  displaying  It  quickly  In  graphical 
or  other  useful  forms  Is  mandatory. 

The  Introduction  of  higher  order  Interpolation  Into  the 
programs  Is  absolutely  necessary  In  order  to  obtain  high  enough 
accuracy  to  solve  engineering  and  scientifically  useful  problems. 
The  preliminary  conclusion  of  Chapter  1  that  the  method  charac¬ 
teristics  is  potentially  the  most  accurate  of  the  three  methods 
considered  there  has  had  no  doubt  cast  upon  It  by  the  results 
of  this  study.  Quite  to  the  contrary.  It  appears  that  rather 
encouraging  results  were  obtained  despite  the  very  crude  linear 
Interpolations  used  to  obtain  the  results. 

The  utility  of  the  method  of  characteristics  applied  to 
multi-dinenslonal  flow  fields  Is  obvious.  In  the  course  of  trying 
to  decide  what  flow  fields  should  be  calculated  as  typical  example 
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cases,  so  many  came  to  mind  that  the  choice  vras  difficult. 

Many  Interesting  unsteady  hypersonic  flows,  which  were  the  moti¬ 
vation  of  this  study,  can  be  calculated  with  the  method.  The 
vertical  oscillation  of  the  blunt  body  presented  In  the  previous 
section  Is  Just  one  Interesting  example.  Other  types  of  oscil¬ 
lations  or  trajislent  motions  can  easily  be  handled  with  the  method. 
Interesting  steady  flows  might  be  determined  by  starting  from 
some  known  steady  flow  and  passing  through  a  specified  transient 
to  a  second  more  Interesting  and  perhaps  previously  unknown 
steady  flow.  The  symmetric  and  asymmetric  warping  results  pre¬ 
sented  In  the  previous  section  give  encouraging  Indications  that 
this  technique  Is  feasible  and  practical.  Many  other  Interest¬ 
ing  applications  of  the  method  can  be  easily  thought  of  by 
merely  letting  the  mind  wander. 

One  limiting  aspect  of  the  method  Is  the  large  amount  of 
Initial  data  required  to  start  the  solution.  Many  theories  are 
available  to  give  certain  specific  details  of  flow  fields,  but 
the  complete  details  of  entire  flow  fields  are  seldom  determined 
and  hardly  ever  published.  For  this  reason  the  application  of 
the  method  might  be  somewhat  limited,  but  there  Is  the  Intrigu¬ 
ing  possibility  of  specializing  the  general  method  In  certain 
cases  so  that  It  could  determine  Its  own  Initial  data.  For 
example,  the  Initial  conditions  for  a  three-dimensional  unsteady 
flow  might  be  obtained  by  using  a  specialized  form  of  the  general 
method  to  solve  a  three-dimensional  steady  supersonic  flow. 
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CHAPTER  6 


EXTENSIONS  AND  FUTURE  DEVELOPMENTS 


6.1  lamtdiate  Improvementi  to  Exiiting  Progrmt 

The  existing  progreos  listed  in  Appendix  C  ere  fsr  from 
being  optimized  with  respect  to  operating  time  end  storage  re-' 
qnirements.  A  few  improvements  which  can  be  made  to  these  pro- 
grasM ,  are  listed  in  this  section.  This  is  not  an  exhaustive 
list  of  improvements ,  but  includes  those  changes  which  appear 
advantageous  from  results  obtained  through  operation  of  the  pro- 
grasis,  up  to  the  present  tism. 

The  iteration  loops  in  the  subroutines  could  be  improved 
by  using  both  residual  tests  and  convergence  test  criteria  in 
the  form  of  (4.1).  Then,  if  the  convergence  test  value  were  too 
small,  the  residual  tests  could  still  release  the  calculation 
frtM  the  loop.  Conversely,  If  the  convergence  test  values  did 
not  require  that  the  residuals  be  reduced  to  zero  to  eight 
decimal  digits,  some  computation  time  could  be  saved  by  releas¬ 
ing  the  calculation  from  the  loop  as  soon  as  the  convergence  cest 
criteria  were  satisfied.  Of  course,  residual  tests  can  be  applied 
to  only  those  loops  which  are  solving  algebraic  equations,  be¬ 
cause  residuals  do  not  exist  for  the  loops  which  are  integrating 
differential  equations.  As  pointed  out  in  Section  3.2,  the 
inn«r  iteration  loops  in  FLDPT,  BDYPT,  and  SHKPT  can  be  elimina¬ 
ted  by  more  closely  utilizing  the  third  order  Richmond  iteration 
scheme.  This  might  reduce  the  computation  time  required  to 
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d^ttniliM  thm  eoordinattt. 

The  •x«cutlv«  prograa  and  THUHET  can  ba  coablnad  into 
ona  program  with  a  langth  alightly  longar  than  TRUMET ' a .  In 
tha  procaaa  of  writing  this  conbinad  program,  tha  possibility 
of  using  highar  ordar  intarpolation  should  ba  invast igatad. 

Hora  labaling  of  tha  printad  output  could  also  ba  includad  in 
tha  program,  togathar  with  BK>re  options  for  supprassing  or  ob¬ 
taining  tha  larga  quantities  of  output. 

Tha  PATCH  subroutina  can  ba  axtansivaly  improvad.  First, 
tha  unsuccassful  option  for  locating  body  and  shock  points  dis¬ 
cussed  in  Section  4.3  can  ba  ramovad,  and  tha  sisa  of  tha  pro¬ 
gram  substantially  reduced.  Tha  iteration  loop  tests  can  ba 
improvad  as  discussad  In  tha  sacond  paragraph  of  this  section. 
Highar  ordar  intarpolation  could  ba  substituted  for  tha  linear 
intarpolation  currently  used  in  tha  subroutina. 

Tha  convergence  criterion  for  tha  components  of  tha  shock 
wave  normal  vector  could  ba  studied  and  perhaps  improvad.  In 
tha  currant  SHRPT  subroutina,  tha  convergence  test  value  has 
bean  arbitrarily  sat  at  tan  times  tha  outer  loop  convergence 
test  value. 

Finally,  tha  executive  and  also,  the  TRWET  programs  could 
ba  much  more  efficiently  written  in  the  FAP  computer  language 
because  they  are  basically  manipulative,  rather  than  arithmetic 
programs.  Of  course,  probably  all  of  the  programs  could  be  im¬ 
proved  somewhat  if  written  by  skilled  programmers  in  FAP. 


6.2  fwtqra  Ptvlop— ntt 

As  nsntionsd  in  ChAptsrs  1  «nd  k,  vhtiisvsr  llBltstiont 
wsrs  iacroduesd  into  ths  problsa  bsing  eonsldsrsd,  thsy  iisrs  dons 
in  such  s  wsp  thst  thsy  could'  bs  aost  ssslly  rsaovsd  st  s  Istsr 
tias.  For  instsncs,  ths  nonorthogonsl  Thornhill  coordlnsts  trsns- 
fdcastion  discussed  in  Section  3*2  could  hsve  been  used  to  sla- 
plify  the  tao-diaensionsl  unstesdy  flow  cslculeted  here.  But, 
it  aes  not  used  beceuse  no  generelixed  trsnsforastion  of  this 
type  is  svsilsble  to  siaplify  the  three-diaensionsl  unsteady 
problea.  In  this  section,  developaents  are  discussed  which  can 
be  added  to  the  prograas  in  the  future.  It  is  felt  that  no 
fuadaaental  changes  will  be  necessary  in  the  prograas,  or  the 
analysis  to  carry  out  these  extensions. 

The  prograas  can  be  extended  to  consider  three -diaens ions 1 
unsteady  flows  when  the  newest  coaputers ,  which  are  being  built 
at  the  present  tiae,  are  available.  These  aachines  will  have 
aagnetic  core  storage  units  of  130,000  to  260,000  words,  and 
will  aaka  the  solution  of  aulti-diaensional  problems  practical. 

The  addition  of  a  single  space  diaension  to  the  prograas  pre* 
seated  here  should  be  rather  siaply  accoaplished.  The  double 
iteration  process  to  solve  the  four  Hach  hyper cone  equations, 
(3.%a),  was  progreaaed  and  debugged  to  check  the  feasibility 
of  the  process  for  four  independent  variables.  The  prograa 
was  written  in  such  a  way  that  the  equations  for  three  Kach 
cones  ia  two-diaensional  unsteady  flow  could  also  be  solved. 

The  prograa  worked  very  well,  but  it  was  concluded  that  its 
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length  urns  Increeted  to  e  lerge  ostent  beceueo  both  throo  end 
four  indepondont  verieblt  eolutione  could  bo  dotoniinod,  rothor 
than  Juat  ona  or  tba  other. 

Changing  the  prograna  to  eonalder  equilibriuB,  real  gaaea 
la  a  rather  alaple  proceie,  as  sMntioned  in  Chapters  3  and  4. 

The  only  practical  drawback  night  be  prohibitive  Increases  In  con- 
putatlon  tlsM.  This  too,  could  be  overcome  by  the  new  computers 
to  be  available  shortly  which  will  use  multiple  processing  methods 
and  equipment  to  decrease  the  time  required  to  complete  a  given 
calculation. 

Monequl librium  flows  can  also  be  attempted  with  these 
higher  speed  machines.  The  Inclusion  of  the  Integration  of  the 
species  continuity  equation  along  the  streamline,  will  require 
major  modifications  to  the  existing  programs,  but  experience  being 
gained  at  the  present  time  in  calculating  two-dimensional  steady 
nonequilibrium  flows,  can  be  utilized  to  facilitate  the  modifica¬ 
tion.  The  increased  ttcrage  capacity  and  higher  speed  of  future 
machines  should  make  the  calculation  of  this  most  general  flow 
possible. 

In  considering  more  general  problems,  points  falling  on 
interior  shock  waves,  contact  surfaces  and  expansion  waves  must 
also  be  calculated.  These  might  result  from  bodies  idilch  have 
comers  and  other  discontinuities.  Procedures  can  be  developed 
to  calculate  such  points  (see  Butler  ,  for  example).  Subroutines 
similar  to  the  FLDPT,  BDYPT,  and  SHKFT  will  have  to  be  written 
for  these  more  complicated  situations.  These  subroutines  appear 
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At  iMtC  fMtlblA,  hVLt  will  rA^lrO  AS  AStADAiVA  prOgTAHiillg 
Affort  to  writA  A  prACtlcAl  tot  of  tatooiitliiAe . 

In  tbo  Boro  dlotAst  futwrii,  it  !•  poAfiblo  thAt  a  unlvoroAl 
lovlicld  flow  progroB  could  bo  dduolopod  which  tolvod  the  threo- 
dlBMitlonAl,  unitAAdy,  lnviicld»  nonAqollibrluB  flow  of  ab  Arblt- 
rAry^Hitsturo  of  gAset  About  a  goner a1  body.  The  progroB  could  bo 
frltton  in  fuch  o  woy  thot  it  could  olto  do  probloBO  which  oro 
•pociol  tinpliflod  cotoc  }of  tho  Boot  gonorol  flew.  At  the  prosont 
tlBo,  thoro  oro  Bony  roodblocko  to  writing  ouch  o  progroB  in  o 
procticol  foTB,  but  thoro  oppoor  to  bo  no  fundoaontAl  rootont  why 
ouch  A  progroB  could  not  bo  crootod. 

It  it  pootiblo  thot  ooBo'dAy,  this  progrAB  could  bo  cob- 
blood  with  other  progroat  thot  conputo'  boundary  layort  or  other 
portlono  of  tho  flow  whore  vlocotlty,  hoot  conduction,  end  dlf- 
futlon  night  bo  laportont.  EloctroBOgnotlc  offoett  night  olto  bo 
inoludod.  At  proiont,  tone  people  ore  thinking  About  utlng  con- 
putort  At  dotlgn  toolt,  to  that  dotlgnort  or  Invottlgatort  working 
At  now  and  yet  to  bo  dovolopod  Input-output  dovicot,  could  Interact 
Boro  directly  with  conputoro  In  tho  dotlgn  procott.  Then,  for 
oxonplo,  on  Aircraft  or  tpacocraft  dotlgnor  could  uto  tho  gonorol 
flow  field  progroBt  dltcuttod  hero,  together  with  ttructurol  ttrott 
onolytlt  progroBt,  hoot  trontfor  progrtat,  control  tyttoat  pro- 
groat  ,  etc.  to  work  ot  a  conputor  cobooIo  end  quickly  end  offl- 
cioatly  dotlgn  tho  vehicle.  Tho  progrtat  would  bo  written  luch 
that  tho  output  of  one  would  bo  utod  at  Input  for  tovorol  othort 
and  vice  vorta.  Thnt,  they  could  bo  utod  In  on  Itoratlvo  aonnor 
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to  optlalso  tho  dofign.  Aftor  tho  vohiolo  it  dttlgntd,  it  it 
ooaetivtblt  that  it  cooli  bo  built  and  aataablad  by  aatoaatad 
■achittt  toola  iihieh  art  diraetly  coimaetad  to  tha  eoaputar  uaad 
in  tha  datign  proeatt. 

It  appaart  that  thara  ara  vast  pottibilitiat  for  tha 
utilisation  of  high  apaad  alactronic  digital  eoaputart  and 
nunarical  solution  nathods  in  solving  coaplieatad  nonlinaar 
problaas.  This  is  aspacially  trua  uhan  vary  ganaral  boundary 
conditions  and  conplic^tad  physical  procassas  nust  ba  conaidarad. 
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APPODIX  A 


iy  ■  - 


Tbm  MthMitlMl  tlMorj  of  ehoroetoristlef  for  hyporbollc 
portiol  dlfforoBtlal  o^Mtiont  is  noil  dovolopod  ond  it  eon- 
tldortd  elottlctl.  Tbo  con^loto  tSMory  will  not  bo  roitorotod 
boro.  0nl7  tboto  proportlot  and  o^notlona  ntilitod  in  dovolop- 
ifl|^  tbo  ^ortienlar  niaorieol  aothod  of  eoneom  to  tbit  tnolytit 
aro  praaantod.  For  tbo  eoaploto  ■atboaatieal  tboory,  tbo  roador 
it  roforrod  to  Chaptor  6  of  tbo  tocond  voluao  of  Conrant  and 
lilbort^^.  Tbo  thoory  appliod  to  tbo  flow  of  a  eoaprottiblo 
flaid  it  ditcuttod  by  von  Mitot^^  (pp.  100-134). 

Ccmtidor  a  tot  of  ^piaai-linoar  firtt-ordor  partial  dif- 
forontial  tifuationt  nbieh  can  bo  mritton  in  tbo  fon* 

A.;,,  +  Bj  -  O  (A.l) 

iX; 

Qoati-linoar  Mont  that  tbo  bigbott-ordor  dorivativot  appoar 
4to  only  tbo  firtt  dogroo.  Koto  that  eontidoring  firtt-ordor 
ognationt  it  not  a  rottriction  bocanto  highor-ordor  tott  of 
agnationt  can  bo  rodncod  to  firtt-ordor  by  dofining  nov 
variablot.  In  (A.l),  i4|^  it  a  gonoral  dopondont  variable,  and 
an  indopondont  variable.  A^j|^  and  Bj  aro  functiont  of 


*  The  index  ttnBttion  notation  of  Cartotian  tontort  it 
onployod. 
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Md  Xi  only.  Kwt,  (A.l)  it  trtnsfoTMd  to  a  mv  tat  of  co- 
ordiaataa, 

it  ° 

Confidor  a  nanlfold  of  on#  last  diaantion  than  tha  parant 
■paca,  a.g.,  for  thraa  indapandant  varlablat,  a  two-dlaans tonal 
aurfaca.  Lat  tha  aanlfold  ba  givan  by 


*  conatant 


(A. 3) 


upon  which  and  all  6  Uy^/  6  axcapt  ^ 

Tha  ^  ^X  ^  dataralnad  by  uaing  aquationa  (A. 2). 

If  (A. 2)  ara  conaidarad  a  conplata  aat  of  algabraic  aquationa 
in  tha  bUy^/  thay  can  ba  aolvad  for  ^(«<|(/  ^  ^X 


Tha  aanifolda  for  which  (A. 4)  ia  not  trua  ara  tha  charactariatic 
■anifolda,  a.g.,  for  thraa  indapandant  variablaa  —  a  charactar¬ 
iatic  aurfaca .  Thua 


(A.5) 


ia  tarvad  tha  charactariatic  aquation  and  dafinaa  tha  charactar¬ 
iatic  nanifoldi.  Mota  that  whan  Uy^  and  ^  ^ 

ara  givan  on  tha  charactariatic  manifold,  b  M^/  ^3^  cannot  ba 
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id«c«niliitd  frcMi  (A. 2)  b«e«u««  of  (A. 5).  This  is  ono  proporty 
of  tho  charactoriitlc  Molfoldi.  Thoro  oro  othor  proportloo 
ifhielr  could  bo  ond  hovo  boon  nood  to  dofino  thoa,  but  thooo  oil 
lood  to  tho  ioao  choroetorittlc  o^uotion,  (A.  3). 

Voxt,  conoldor  o  llnoor  eoablnotlon  of  tho  oquotiono  (A. 2), 
uhleh  COB  bo  vritton  in  tho  fora 


>j  Aijk  ^ 


or  thlo  con  oloo  bo  vrltton  in  tho  fora 


(A.6) 


Ck  +  0^ 

Vi 


uhoro  tho  doflnitlon  of  tho  now  dyabo'io  io 


,+  &-0  (A.7) 


‘I.  ■  'A 

(A.8«) 

\  -  -li  j(l\ 

(A.  8b) 

H  -  A; 

(A. 8c) 

\  -  Aj  Aijk 

• 

• 

(A.8<i) 

• 

• 

• 

c  -  Aj  Bj 

• 

(A.8«) 

120 


On  thn  eharnetnrittic  annlfold  glvnn  by 


*  eonstnnt 

(A. 3)  holds,  so  Aj  esn  bs  found  such  thst 


=  Aj  Aiji^ 


(A.9) 


Thsn,  on  ths  chsrsctsrlstic  Mnifold,  oqtistion  (A.  7)  has  tbs 
font 


-I-  fi  =  O  (A.IO) 


%rh«rs  D^,  E^,  and  G  are  given  by  (A. 8b)  through  (A. 8a)  and 
A  j  is  determined  from  (A.9).  Equation  (A.IO)  is  called  the 
coapatibility  equation  and  has  the  useful  property  that  deriva¬ 
tives  normal  to  the  characteristic  manifold  do  not  appear.  In 
general,  there  might  be  more  than  one  compatibility  equation  for 
a  particular  characteristic  manifold,  but  there  must  be  at  least 
one  such  equation.  The  compatibility  equations  are  basic  for  the 
method  of  characteristics. 
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APPENDIX  B 


A  Note  on  Stability 

Consider  a  set  of  first  order  hyperbolic  pertiel  dif 
ferentisl  equetions  of  the  form 


mtiere  f  is  e  vector  mede  up  of  n  dependent  verieblee,  (K^  ere 

■  specs  coordinetee  end  Ar  ere  reel  constent ,  coefficient  metricee 

52 

Covrsnt,  Friedrichs  end  Lewy*^  determined  thet  e  necessery  con¬ 
dition  for  convergence  of  e  difference  scheme  for  this  set  of 
equetions  is  thet  the  domein  of  dependence  of  the  dlffrence 

scheme  must  contein  the  domein  of  dependence  of  the  differentiel 

59  524 

eqnetion.  Hehn*^’’  ,  using  the  work  of  Lex'^  ,  showed  thet  for 
simplicisi  difference  schemes  (i.e. ,  schemes  thet  use  e  mlnimimn 
nwber  of  points  in  the  initiel  surf  see  to  determine  e  new  point) 
the  Coursnt-FriedriciS'lmwy condition  is  sufficient  es  well  es 
necessery  for  convergence,  end  therefore,  elso  for  stebility. 

The  lineerised  equetions  for  inviscid  nonequilibrium 
thrme-diswnsionel  unsteedy  flow  cen  be  put  in  the  form  of  eque- 
tion  (1.1).  The  equetion  for  one-  end  two-dimensionel  unsteedy 
flow  end  two-  end  three -dimensionel  supersonic  steedy  flow  elso 
heve  this  form.  Hence,  the  stebility  end  convergence  criterion 
mentioned  in  the  ebove  peregreph  epplles,  et  leest  locelly 
(becsttse  of  the  lineerlsetion) ,  to  these  types  of  flows. 
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It  it  poffiblc  to  rtmovc  the  restriction  thst  the  criterion 


BS. 


applies  only  locally  because  of  the  linearisation  if  the  equa* 

tions  and  the  true  solution  possess  a  certain  degree  of  saoothness. 

55 

Strang  has  sho%m  that,  if  the  equations  and  their  solution  possess 
enough  continuous  derivatives,  then  the  convergence  of  the  full 
nonlinear  equations  depends  on  the  stability  of  the  linearized 
equations.  The  existence  of  the  continuous  derivatives  of  the 
solution  cannot  be  proved  before  the  solution  is  known,  but  it 
might  be  expected  that  for  "reasonably  well  behaved"  physical 
problems,  the  solution  will  possess  the  required  degree  of 
smoothness . 

This  stability  criterion  can  be  applied  to  the  basic 
finite  difference  networks  proposed  for  problems  with  three  in¬ 
dependent  variables,  and  discussed  in  Section  3>1*  The  tetra¬ 
hedral  characteristic  line  network  of  Figure  3  is  found  to  be 
unstable  trhen  the  stability  criterion  is  applied.  A  good  approx¬ 
imation  to  the  domain  of  dependence  of  point  in  the  initial 
surface,  is  a  circle  through  points  P2 ,  and  P^.  This  is  the 
domain  of  dependence  for  the  partial  differential  eqiutions. 

The  domain  of  dependence  of  the  difference  net  is  the  triangle 
with  vertices  at  P^ ,  P2 ,  and  P^  (formally  termed  the  convex  hull 
of  the  difference  scheme).  The  domain  of  dependence  of  the  dif¬ 
ference  scheme  does  not  contain  the  domain  of  dependence  of  the 
differential  equations,  and  therefore,  does  not  fulfill  the 
stability  criterion. 
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PhyticAlly,  Ch«  acability  criterion  teemt  very  logical, 
becauae  it  ia  not  hard  to  imagine  that  in  certain  areaa  of  a 
flow  field,  pointa  lying  within  the  circle  but  outaide  the  tri¬ 
angle  circumacribed  by  the  circle ,  could  have  a  large  influence 
on  the  new  point  being  calculated.  However,  in  thia  network, 
effeeta  of  pointa  outaide  the  triangle  are  not  taken  into  account 
in  the  numerical  calculation,  ao  that  the  aolution  could  be 
erroneoua . 

To  inaure  atability,  a  modification  to  the  tetrahedral 
diaracteriatic  line  network  (ahown  in  Figure  6)  ia  propoaed. 
Starting  with  the  original  baae  pointa  Pj^,  P2 ,  and  P^,  inter¬ 
mediate  baae  pointa  Pj^^  *  ^23*  ^13  determined  by  in- 

acrlblng  a  circle  within  the  triangle  with  verticea  at  the 
baae  pointa.  The  pointa  of  tangency  between  the  circle  and 
the  aldea  of  the  triangle  are  choaen  aa  intermediate  baae  pointa 
and  the  flow  proper tiea  at  theae  pointa  are  determined  by  inter¬ 
polation.  The  intermediate  baae  pointa  are  uaed  in  exactly  the 
aame  way  aa  the  baae  pointa,  P^^,  P2»  and  P^,  were  uaed  in  the 
original  network.  Thia  modified  network  fulfilla  the  atability 
criterion,  in  that  the  circular  domain  of  dependence  of  the  dif¬ 
ferential  equationa  falla  within  the  triangular  domain  of  depend¬ 
ence  of  the  difference  acheme.  If  the  aame  baae  pointa  are  uaed 
to  dalculate  a  new  point  with  both  the  modified  and  unmodified 
networka ,  the  atep  aize  from  the  initial  aurface  to  the  new 
point  ia  amaller  for  the  modified  network.  Thua,  the  uae  of  the 
modified  network  requirea  a  larger  nwber  of  calculationa  to  do 
a  given  problem,  but  it  ia  well  known  that,  in  atandard  finite 


124 


dlffor«nc«  techniques,  decreesing  the  step  size  has  e  stebi- 
lising  influence. 

A  niSMricel  celculetion  has  been  nuide  of  the  flow  of 
a  perfect  gas  about  a  circular  cylinder  with  Its  axis  normal 
to  a  Mach  3  free  stream.  Both  the  modified  and  unmodified  net* 
works  were  used  to  calculate  the  flow  with  the  body  surface  held 
constant  in  time.  It  was  found  that  the  pressure  along  the  axis 
of  symmetry  was  the  most  sensitive  indicator  of  the  onset  of 
instability.  Figure  33  shows  plots  of  the  pressure  ratio  versus 
the  coordinate  parallel  to  the  free  stream  velocity  for  points 
closest  to  the  axis  after  various  time  steps  in  the  calculation. 
Network  A  is  the  unmodified  tetrahedral  characteristic  line  net¬ 
work,  while  Network  B  is  the  modified  form  of  the  same  net.  Note 
that  these  points  are  not  located  exactly  on  the  axis  of  symmetry 
and  do  not  have  exactly  the  same  space  coordinates  at  each  step, 
so  that  the  magnitudes  of  the  pressure  ratio  should  not  be  exactly 
the  same  on  each  plot,  but  they  do  very  graphically  indicate  the 
onset  of  instability.  The  instability  with  the  unmodified  network 
had  grown  so  large  by  the  eighth  step  that  calculation  could  not 
be  continued,  whereas  no  instability  has  been  detected  idiile  using 
the  modified  network. 

The  tetrahedral  characteristic  surface  network  is  also 

simplicial  and  is  found  to  be  stable  when  the  stability  criterion 

29 

is  applied.  Tsung  utilized  this  network  and  indeed,  found  no 
instability  in  his  calculations. 
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The  left  three  networks  discussed  in  Chapter  3  and  shown 
in  Figures  8  to  10,  are  not  siaiplicial  and  the  stability  criterion 
cannot  be  foxmally  applied.  The  physical  arguaent  given  in  the 
fifth  paragraph  of  this  Appendix  would  still  seem  to  apply  to 
these  networks,  however.  Hence,  application  of  the  criterion 
would  seen  to  give  some  insight  into  the  stability  of  the  net , 
even  though  it  does  not  formally  apply. 

It  appears  that  these  three  networks  might  be  unstable, 

if  only  the  local  net  configuration  is  considered.  However, 

-JQ 

Moretti,  et.al.-^  has  utilized  the  network  of  intersections  of 

28 

reference  planes  with  characteristic  surfaces,  and  Butler  ,  the 
pentahedral  bicharacteristic  line  network,  and  they  encountered 
no  instability  in  their  calculations.  Closer  examination  of  the 
networks  reveals  that  in  each  of  the  calculations,  the  base  points 
in  the  initial  surface  must  be  moved  about  in  the  initial  surface 
as  the  solution  at  the  new  point  is  obtained  by  iteration.  This 
requires  that  the  properties  at  the  base  points  must  be  obtained 
by  interpolation  in  the  initial  surface.  The  interpolation  schemes 
uaed  by  Moretti  and  Butler  are  such  that  they  effectively  increase 
the  domain  of  dependence  of  the  difference  scheme  to  the  extent 
that  it  contains  the  domain  of  dependence  of  the  differential 
equations.  Hence,  the  overall  schemes  could  be  stable. 

Now,  in  applying  the  stability  criterion  to  the  problem 
of  four  independent  variables,  four  base  points  in  an  initial 
hyparsurface  must  be  considered.  A  net  utilizing  four  base 
points  is  simplicial  so  the  criterion  can  be  applied.  The  four 
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b«s«  points  dsfins  s  tstrshcdron  in  the  initial  hTper surf act 
iifhich  is  the  domain  of  dspandanca  of  the  finite  difference 
schesie.  The  domain  of  dependence  of  the  partial  differential 
equations  is  the  intersection  of  the  characteristic  hTpersurface 
(the  Mach  hyperconoid)  with  the  initial  hypersurface.  If  the 
initial  hypersurface  is  a  t '  -  constant  hypersurface,  the  inter* 
section  is  a  sphere  or  can  be  closely  approximated  with  a  sphere 
for  most  reasonable  initial  hypersurfaces.  Hence,  if  the  domain 
of  dependence  of  the  finite  difference  scheme  is  to  contain  the 
domain  of  dependence  of  the  differential  equations,  the  sphere 
must  be  contained  within  the  tetrahedron.  The  modified  tetra¬ 
hedral  characteristic  line  network,  shown  in  Figure  11  and  dis¬ 
cussed  in  detail  in  Chapter  3,  fulfills  the  above  requirement  and 
therefore,  should  give  stable  results. 

Previously,  it  has  been  tacitly  assisMd  in  characteristic 
calculations  involving  three  independent  variables  that  inte¬ 
gration  on  the  characteristic  surface,  or  at  least  along  bi¬ 
characteristic  lines,  assures  the  convergence  and  stability  of 
the  process.  This  is  not  the  case.  Care  must  be  taken  in 
choosing  a  net  configuration  which  is  stable  and  convergent.  The 
Courant-Friedr ichs-Lewy  condition  can  be  employed  as  a  test  for 
stability  and  convergence  of  a  finite  difference  net%rork  which 
is  simplicial,  and  it  can  be  argued,  physically,  that  this  pro¬ 
vides  a  necessary  condition  for  the  stability  of  networks  %rhich 
are  not  simplicial. 
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APPEMDIX  C 


Progr—  Flow  Dia£rMi»  and  Ligtintt 

It  li  felt  that  ezect  program  listingt  ere  mich  more  uee> 
ful  then  deteiled  flow  diegreae  in  etteaptlng  to  duplicete  or 
extend  exieting  coaputer  progreae.  Deteiled  flow  diegreai  still 
require  auch  time  end  effort  to  trenslete  thea  into  working  coa- 
puter  progreae,  end  beceuee  of  their  deteil,  they  ere  not  eesy 
to  use  in  etteapting  to  understand  the  broed  pattern  in  the 
progrea.  Listings,  on  the  other  hand,  require  only  keypunching 
to  obtain  working  progreas.  For  this  reason,  only  a  few  very 
general  flow  diegreas  for  the  FLDPT,  BOTPT,  end  SHKPT  subroutines 
ere  included  here  to  give  a  general  idee  of  the  progrea  organisa¬ 
tions.  Exact  copies  of  the  listings  are  presented  to  aake  the 
prograas  easily  reproducible.  These  listings  are  not  in  the 
neatest  possible  fora,  but  have  been  left  in  the  exact  fora  to 
which  they  evolved  in  the  debugging  process.  The  listings  have 
been  left  in  this  fora  to  eliainate  the  possibility  of  intro¬ 
ducing  transcriptions  errors  in  a  process  of  putting  thea  in  a 
neater  foraat. 

The  M.I.T.  Coaputation  Center  has  altered  its  systea,  so 
that  READ,  PRINT,  and  PUNCH  FORTRAN  stateaents  read  and  write 
tapes  and  thus  all  input  and  output  is  done  off-line.  XSIMEQF 
is  a  FAP  coded  subprograa  avai labia  on  the  Computation  Center 
Library  tape  which  solves  the  aatrix  eq\iation 

(A)  [X]  -  [B] 
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whcr«  A  and  B  are  given  and  X  la  deterninad.  EXIT  la  Che  aub- 
roucina  called  to  acconpllah  a  nontal  CarBlnacion  of  the  progr 
execution . 

The  flow  dlagrama  and  llatinga  follow  in  the  following 

order : 


Diacrama 


Llatinga 


Flow  Diagram  for  the  FLDPT  and  BDYPT  Subroutinea 
Flow  Diagram  for  the  SHKPT  Subroutine 

Executive  Program 

SCAN 

ADDROW 

TRUNET 

TINTER 

PATCH 

FLDPT 

SHKPT 

BDYPT 

Several  XCHANG'a 

General  form  of  XCHANG 

XINTER 

ADET 

ROEDET 

THETA 


129 


START 


CALCULATE 
FIRST  ESTIMATES 


calculate  the 

COORDINATE  CORRECTIONS 


HAVE  THE 

COORDINATE  CORRECTIONS 
CONVERGED 


CALCULATE 
THE  COORDINATES 


HAVE  ThC^ 
COORDINATES 
CONVERGED^ 


CALCULATE  THE 
FLOW  PROPERTIES 


<^HAVE  THE 
FLOW  PROPERrriES' 
CONVERGED^ 


RETURN 


FLOW  DIAGRAM  FOR  THE  FLDPT 
AND  BDYPT  SUBROUTINES 
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I 


START 


CALCULATE 
FIRST  ESTIMATE 


CALCULATE  THE 
COORDINATE  CORRECTIONS 


HAVE  THE 

COORDINATE  CORRECTIONS 

VE  RGED,,^^''^ 


CALCULATE 
THE  COORDINATES 


HAVE  THE  ^ 
COORDINATES 
CONVERGED 


CALCULATE 
FLOW  PROPERTIES 


CALCULATE  CORRECTION 
TO  THE  SHOCK  WAVE 
NORMAL  VECTOR 


HAVE  The 

FLOW  PROPERTIES  AND 

normal  vector 

CONVERGED 

■7 

ivts 


return) 


FLOW  DIAGRAM  FOR 
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TABLE  I 


Rtiults  of  a  Typical  Fitld  Point  Calculation 


Input  Points 


(*,  y.  t'  in 
BTU/lbo.  OR) 


fps,  p  in  pif, 


lit  point 
2nd  point 
3rd  point 


-1.1447676 

-1.0939953 

-1.1141122 


0.58913501 
0. 54756027 

0.49530905 


0.14216991 

0.14148745 

0.14097986 


1st  point 
2nd  point 
3rd  point 


u 

874.92389 

727.29974 

648.66292 


V 

897.04218 

923.78306 

849.49358 


E 

6071.4579 

6218.8732 

6405.26^3 


s 

3.1160775 

3.1168170 

3.1169245 


PATCH  Innar  Loop  Middlt  Loop  Outar  Loo] 

0.0001  0.0001  0.0000005  0.0001 


Coordinatas  of  Points  Outi 


PATCH 


-1.1018552  0.5545i524 
-1.0981097  0.53687342 
-1.1279440  0.53764377 


t’ 

0.141^654 

0.14138363 

0.14151682 


First  Estimata  of  tha 

Coordinatas 

X 

I 

t;. 

-1.1079344 

0.5487763 

0.1552792 

First  Innar  Loop  Iteration 

Ax 

At' 

-1.6167161xl0‘lt 

4.5400817x10*: 

1.1548015x10*1! 

-1.5992165x10*7 

4.5651031x10*: 

1.2058911x10*7 

-1.5992875x10*: 

4 . 5650098x10*: 

1.2056934x10*7 

-1.5992875x10*^ 

4.5650102x10*^ 

1.2056942x10*^ 
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Firit  Middle  Loop  Raault  (no  iteration  required) 

11  t ' 

-1.1080943  0.5492328  0.1553947 

(Residiialf  all  leas  than  5  x  10*^) 


First  Outer  Loop  Result 


E 


s 


+  V 


2 


726.8993  912.92059  6171.1549  3.1168188 


1166.9651 


Second  Inner  Loop  Iteration 


7.2747^4x10"! 

7.2885107x10"! 

7.2885356x10"^ 


1.2188^92x10 

1.2169935x10 

1.2169902x10 


-4 

-4 

-4 


-5.29048^x10"! 

-5.2511849x10"! 

-5.2511145x10"^ 


Second  Middle  Loop  Result 

-1.1080214  0.545i111 


Second  Outer  Loop  Result 

u  y  £ 

725.35434  915. 54065  6175.2055 


Third  Inner  Loop  Iteration 

^  ^  .A 

-5.2475364x10  I  8.5396863x10  ” 
-5.2465308x10"®  8.5410568x10"® 

-5.2465308x10"®  6 . 5410568x10'® 


t ' 

0.15^422 


8  \/u^~T”v^ 

3.1188192  1168. O555 


-1.3029714x10"? 

-1.3001144x10"? 

-1.3001143x10"® 


Third  Midd le  Loop  Result 

-1.1080266  0.54^1196 


t ' 

0.155^409 


Third  Outer  Loop  Result 
u  V 


725.26143  915. 86852 


E  • 

6176.1796  3.1188191 


\/u^  + 
1188.25^8 
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Fotxrth  Inner  Loop  Iteration 


-6.2957^^56x10’? 

-6.2954510x10’; 

-6.2954509x10’' 


1.59815578x10’? 

1.5981327x10’? 

1.5981327x10"® 


2.469?757x10’? 

2.4702034x10’; 

2.4702034x10’' 


Foiirth  Middle  Loop  Remit 


-1.10F0272 


0.545i2119 


t ' 

0.1553^113 


Fourth  Outer  Loop  Rciulf 
u  V 

725-24378  915-91642 


£ 

6176.2986 


3.1168191 


/u^  + 

1168.2814 


Th«  SABt  point  was  re-calculated  with  the  following  Convergence 
Test  Valuea. 

PATCH  Inner  Loop  Middle  Loop  Outer  Loop 

0.00001  0.00001  0.0000001  0.000001 

The  coordinates  fron  PATCH  and  the  first  four  iterations  through 
the  outer  loop  were  exactly  the  same  as  given  above.  Due  to  the 
decrease  of  the  outer  loop  convergence  test  value  by  a  factor  of 
100,  the  outer  loop  was  iterated  two  more  cycles  with  the  follow¬ 
ing  results. 


Fifth  Outer  Loop  Results 

-1.1080286  0.5495215 

U  V  £ 

725-24043  915-92667  6176.3179 


t ' 

0.l55Pr07 

s 

3.1108191 


iT5F:IB74 


Sixth  Outer  Loop  Results 

-1.1000277  0.54952081 

U  V  £ 

725-73881  915-02730  6176.3199 


t' 

0.15574049 

8  ^u^  +  v^ 

3.1108191  1T0072008 
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TABLE  II 


TIME  COORDINATES  OP  THE  TIME  SURFACES 
FOR  THE  POUR  BLUNT  BODY  CASES 


t  -  t'/U  ■  t'/2000,  t  In  seconds  and  t'  In  feet 

An  asterisk  Indicates  a  time  surface  on  which  ADDROW 
was  used  to  add  points  by  Interpolation 


Time 

Surface 

Number 

Steady 

Flow 
t ' 

Symmetric 

Warping 

t' 

Asymmetric 
Warping 
t ' 

Vertical 

Oscillation 

t' 

0 

.0 

.0 

.0 

.0 

2 

.0858  -  .0182 

.0858  -  .0179 

.0858  -  .0179 

.0858  -  .0182 

4 

.0455 

.0426 

.0431 

.04520 

6 

.1268  -  .0645 

.1237  -  .0603 

,1239  -  .0609 

.1273  -  .0642 

8 

.0930* 

.0838* 

0857 

.0912* 

10 

.1650  -  .1076 

.1555  -  .0980 

.1568  -  .1001 

.1652  -  .1059 

12 

.1339 

.1202 

.1227 

.1327 

14 

[.2033  -  .1493 

.1884  -  .1355 

.1899  -  .1381 

.2045  -  .1483 

16 

.1756* 

.1606* 

.1637* 

.1744 

18 

.2406  -  .1888 

.2235  -  .1744 

.2253  -  .1774 

.2425  -  .1905 

20 

.2153 

.2001 

.2028 

.2174* 

22 

.2372 

1 

.2606  -  .2145 

.2622  -  .2172 

.2802  -  .2309 

24 

.2980  -  .2513 

.2387* 

.2415* 

.2569 

26 

1  .2773* 

1 

.2956  -  .2507 

.2972  -  .2539 

.3178  -  .2709 

28 

.3364  -  .2896 

.2748 

.3519  -  .2689 

.2968* 

30 

.3156 

.3301  -  .  2872 

.2955* 

.3554  -  .3085 

32 

.3735  -  .3284 

.3108 

.3489  -  .3065 

.3359 

3^ 

.3542 

.3647  -  .3236 

.3303 

.3940  -  .3483 
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TABLE  II  (continued) 


. 3759* 

.4304  -  .3867 

.4124 

.4341 

.4556* 

.5091  -  .4648 
.4904 

.5487  -  .5000 


.3469* 

•3986  -  .3573 

.3808 

.4318  -  .3915 

.4146 

.4646  -  .4255 
.4482 


.3824  -  .3416 
.3650 

.4165  -  .3767 

.3997 

.4504  -  .4116 

.4344* 

.4824  -  .4437 
.4663 


.3748 

.4321  -  .3880 
.4140* 

.3690  -  .4243 

.4494 

.5035  -  .4601 
.4853 

.5390  -  .4960 

.5214 

.5758  -  .5321 
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Figur*  1  Charocttrittic  ntt  for  two  indopondtnt  vorioble& 


bichorocftrittic  line 
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1 


domain  of 
d*r«n<J#nc#  of  P4 

OomO'n  of  d«p«ndtnct 
of  tn«  d>ff«rtnct 
tcn«m« 


Figure  5  Tetrahedral  Characteristic  Line  Network 


Figure  6  Modified  Tetrihedral  Charnc terlstlc  Line  Network 


*  * 


Fifurt  6  Nttvorh  of  intorooctioo  of  rofortnco  plono«  wtm  chorocttnstic  luffocot 


X60 


lint  nttwtrk  for  four 


Flfuroll  Wkprontnfotiin  of  tfit  l.lr.ho^ 

iwiiiowiont  votIoMm. 


r  3 


Ib^ 


4 

k 


•hoek--^^ 

hyptrturftct 


hyp«rplon«  normol  to  N4 
powing  through  P4 


Figuro  19  floprttontotion  of  0  portion  of  tho  ohock  point  notwork. 


Figuro  16  RoproMototion  of  0  mothod  for  ditormlning  ttoody  flout. 
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ihocli  hyptrturfoc* 


body  hyporourfoco 


Figurt  17  RtprotoMotlon  of  •  toeoiid  mottiod  for  dolormininf  ttoody  f  o««. 


I 


ttffCulivt  I  inttrpototton  ond  tfoto  twbfoulmtt  which  I  ounliory  mttr 

program  hondlmg  tobroulinci  coicuioU  n««  doto  polotion  lutfOvitmt 


TigurtlS  SatKOu'm*  coiling  ttqufnct 


#«•)!  turfact  Oawtat  aM  aaiatafai  tarlaia  tiiaaw 

arftwtaitaN  la  catcaiaia  aaat  arfaaiaatiw  la  caiciiiaia  aaii 


Tha  iMaa  piMt  ai  laa  aaritcat  al  aaca  al  iha  triaafiaa  ara  aaai  la  eaUaiaia  a  t 
Hmi  aa  iha  aaii  aarfaea 


Ftfvra  19  Caatuiiaa  pra9>aai  a»taa>»aHaa 


aa*ai  laiarpalaiiaa 


Fiawrt  20  twia»palatiaa  a( 


Ma9  9f  MTCM  ia  Ilia  Mtiai  aarfaca 


U0«  4864.1  fpt 
Vco  *  0 

PcD  *  242.2  ptf 
tcD  *  2.926  BTU/lbm.  *R 
y  ■  1.4 

FOR  A  PERFECT  GAS 
1/2 »4240  psf 
ALTITUDE  SOjOOO  ft. 


SHOCK  WAVE 


SONIC  LINE 


o  +  \ 
o  +  a 

^  o  -*-1 
►  o  -»• 

0  4-0 


I  ^  ^  °  ^  ® 

Mr.  ®  ^  Jt 


BODY  SURFACE 


o  BELOTSERKOVSKII 
DATA  POINTS 


X  POINTS  ADDED  WITH 
LINEAR 

INTERPOLATION 


-1.6 


fir  o  if  ^  i  LINEAR 

,  "  O  »  /  INTERPOLATION  / 

*oxoxox<4  / 

>-x^o-x— o-^x-o-k-<4— — J ^  1  — 

-14  -1.2  -1.0  -0.0  -0.6  -04  -0.2 


I 


FIGURE  21  BLUNT  BODY  INITIAL  DATA  NET 
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FIGURE  22  BODY  SURFACE  PRESSURE  FOR  STEADY  FLOW 
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FIGURE  23  BODY  SURFACE  u  VELOCITY  COMPONENT  FOR  STEADY  FLOW 


-030 
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FIGURE  24  BODY  SURFACE  v  VELOCITY  COMPONENT  FOR  STEADY  FLOW 
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FIGURE  25  PRESSURE  ON  THE  MIDDLE  Rl 


FIGURE  26  PRESSURE  ALONG  THE  AXIS  OF  SYMMETRY  FOR 
STEADY  FLOW 


2 
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FIGURE  27  BODY  SURFACE  PRESSURE  FOR  SYMMETRIC  WARPING 


x/R 

FIGURE  28  (o) 

P 


-15  -14  -1.3  -1.2  -II  -1.0 

x/R 

FIGURE  28  (b» 


FIGURE  28  PRESSURE  ON  THE  AXIS  OF 
SYMMETRY  FOR  SYMMETRIC 
WARPING 
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,/R 

FIGURE  29(6) 

1_  f  ' 


t.h 

FIGURE  29(c) 

FIGURE  29  BODY  SURFACE  PRESSURE  FOR  ASYMMETRIC  WARPING 
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BODY  MOTION  ESSENTIALLY  STOPPED, 


(0) 


(b) 


FIGURE  30  VELOCITY  FIELDS  FOR  ASYMMETRIC  WARPING 


176 


figure  31(a)  BODY  SURFACE  PRESSURE  FOR  THE 


E  31  (c  )  BODY 


*2^“® 


Ido 


FIGURE  31  (d)  BODY  SURFACE  PRESSURE  FOR  THE  VERTICAL  OSCILLATION 


2  Ao“co 
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FIGURE  31(e)  BODY  SURFACE  PRESSURE  FOR  THE  VERTICAL  OSCILLATION 
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FIGURf.  32.  TYPICAL  INITIAL  SURFACE  NET  CONFIGURATION  FOR  WEDGE  FLOW 


p  IN  psf  FROM  OBLIQUE 
SHOCK  EQUATIONS 


20,000  h 


26,00  0  h 


PERFECT 

AGREEMENT 


24,000  h 


BODY  SURFACE 
DISCONTINUITY 
TOO  LARGE 


24,000  26,000  28,000 

p  IN  psf  FROM  THE  METHOD  OF  CHARACTERISTICS 

(o) 


WEDGE  ANGLE  CHANGE 
(b) 

FIGURE  33  WEDGE  FLOW  RESULTS 


18^ 


24.8 
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FIGURE  34  TYPICAL  PRESSURE  DISTRIBUTIONS  ON  THE  WEDGE  SURFACE 


‘■p/Poo 

-  INITIAL  CONDITIONS 


AFTER  SIX  TIME  STEPS 


Figure  33  Pressure  Near  the  Axis  of  Synmetry  for  the 
Blunt  Body  Steady  Flow 
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